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Abstract 

We calculate the neutral kaon mixing parameter Bk in unquenched lattice QCD using asqtad- 
improved staggered sea quarks and domain- wall valence quarks. We use the "2+1" flavor gauge 
configurations generated by the MILC Collaboration, and simulate with multiple valence and sea 
quark masses at two lattice spacings of a ~ 0.12 fm and a ~ 0.09 fm. We match the lattice determi- 
nation of Bk to the continuum value using the nonperturbative method of Rome-Southampton, and 
extrapolate Bk to the continuum and physical quark masses using mixed action chiral perturbation 
theory. The "mixed-action" method enables us to control all sources of systematic uncertainty and 
therefore to precisely determine Bk', we find a value of B'P''^'''^{2GeV) = 0.527(6)(21), where the 
first error is statistical and the second is systematic. 
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I. INTRODUCTION 



The kaon B-parameter (Bk), which parametrizes the hadronic part of CP- violating neu- 
tral kaon mixing, plays an important role in flavor-physics phenomenology. When combined 
with an experimental measurement of indirect CP- violation in the kaon sector, ex, Bk con- 
strains the apex of the Cabibbo-Kobayashi-Maskawa (CKM) unitarity triangle. Because 
is known to sub-percent accuracy, this constraint is limited by the theoretical uncertainties 
in several quantities, including Bk- Physics beyond the standard model generically predicts 
additional quark flavor-changing interactions and CP- violating phases. These will manifest 
themselves as apparent inconsistencies between measurements that are predicted to be iden- 
tical within the framework of the standard model. Thus precise experimental measurements 
of quark-fiavor changing weak-interaction processes are sensitive probes of new physics, pro- 
vided that the corresponding theoretical calculations are also sufficiently precise. In this 
work we calculate B^ using lattice QCD with all sources of systematic uncertainty under 
control. This result is needed to interpret the experimental measurement of ex as a con- 
straint on the CKM unitarity triangle, and hence to constrain physics beyond the standard 
model. 

Because an accurate determination of Bk is essential for flavor-physics phenomenology, 
many lattice QCD calculations of Bk have been done over the past decade, each improving 
upon the previous one. The benchmark calculation by the JLQCD Collaboration contains 
a thorough study of the quark mass and lattice spacing dependence [3]. Because it does 
not include the effect of sea quark loops, however, the final result for Bk has a quenching 
uncertainty which is difficult to estimate. The HPQCD Collaboration eliminates this source 
of error in Bk by using dynamical staggered fermions at a single lattice spacing The 
additional species of staggered quarks, referred to as "tastes", however, complicate the 
lattice-to-continuum operator matching procedure, and lead to a ~ 20% systematic error in 
Bk due to neglected higher-order operators and mixing with operators specific to staggered 
fermions that break fiavor symmetry. The RBC and UKQCD Collaborations' calculation 
of Bk contains the effects of three flavors of dynamical domain-wall fermions and employs 
nonperturbative operator renormalization j^. Although their result for Bk has a ~ 6% total 
uncertainty, it relies on a single lattice spacing and an estimate of the size of discretization 
errors based on the earlier quenched calculation in Ref. 4]- 
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Our mixed- act ion lattice QCD calculation combines domain-wall valence quarks and stag- 
gered sea quarks, following the method of the LHP Collaboration We use the "2+1" fla- 
vor asqtad-improved staggered lattices generated by the MILC Collaboration, which include 
the effects of three light dynamical quarks js]. These configurations are publicly available 
with a large range of quark masses, lattice spacings, and volumes and allow for good control 
over the systematic error from chiral and continuum extrapolation [3]. We generate domain- 
wall valence quark propagators using the Chroma lattice QCD software package 8|. The 
approximate chiral symmetry of domain-wall quarks simplifies both the extrapolations to 
physical quark masses and zero lattice spacing and the lattice-to-continuum operator match- 
ing. Because the mixed-action AS* = 2 lattice operator used to calculate the matrix 
element is composed of domain-wall valence quarks, which do not carry the taste quantum 
number, it only mixes with other operators of wrong chirality (due to small residual chiral 
symmetry breaking), not incorrect tastes. This makes the relevant mixed action chiral per- 
turbation theory (MAxPT) more continuum-like than in the purely staggered case. There 
is only one new parameter in the one-loop MA^FT expression for Bk with respect to the 
purely domain- wall case, and it can easily be determined from the staggered pseudoscalar 
meson mass spectrum [3, Use of domain- wall valence quarks in the AS = 2 lattice 
operator also makes the nonperturbative operator matching via the Rome- Southampton 
method [l^ as simple as in the purely domain-wall case. Thus the mixed action method 
combines the advantages of staggered and domain-wall fermions without suffering from their 
primary disadvantages and is well-suited to the calculation of Bk- 

The MILC gauge configurations make use of the rooting procedure to remove the addi- 
tional staggered quark species (tastes) from the calculation of the fermion determinant. A - 
though the "fourth- root trick" has not been proven correct, both theoretical arguments fll- 



and numerical simulations 



15 



17| support the validity of the rooting procedure. Most 



of this evidence is summarized in reviews by Diirr, Sharpe, Kronfeld, and Golterman [18|- 
I21I 1. Given the wealth of evidence substantiating the fourth-root trick, we work under the 
plausible assumption that the continuum limit of the rooted staggered theory is QCD. 

Our calculation of B^ relies upon the ability to correctly extrapolate to the physical quark 
masses and zero lattice spacing using MA^PT 22|, which describes the pseudo-Goldstone 
boson sector of the mixed-action lattice theory. In Ref. [l7| we have therefore performed 
a strong check of the ability of MAXPT to accurately describe the quark-mass and lattice- 
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spacing dependence of the isovector scalar correlator. The ao correlator is particularly 
sensitive to unitarity-violating discretization effects in the mixed-action theory because it 
receives contributions from flavor-neutral two-meson intermediate states. At next-to-leading 
order (NLO), the size and shape of these "bubble" contributions to the scalar correlator are 

n 

completely predicted within MA^PT [15|, given knowledge of a few low-energy constants 
that are easily determined in fits to pseudoscalar meson mass data. We find that, for all 
valence-sea mass combinations on both the coarse and fine lattices, the MAXPT prediction 
is in good quantitative agreement with the numerical lattice data, despite the numerically 
large discretization effects due to the staggered sea sector. Thus we conclude that MAXPT 
describes the dominant unitarity-violating effects in mixed-action lattice simulations. For 
the case of most weak-matrix elements, including B^, NLO MA^PT predicts that unitarity- 
violating discretization effects in 1-loop chiral logarithms are below a percent on the coarse 
and fine MILC lattices 9|. This fact, in conjunction with our successful analysis of the 
scalar correlator, substantiates the claim that we can use MAXPT to remove these effects 
from Bk and to precisely determine its value in the continuum. 

We have also performed a more general check of our ability to control systematic errors in 



jy calculating the light pseudoscalar meson decay 



our mixed-action numerical simulations 

constants, /^r and fx, and their ratio |23|. We use the same gauge configurations and 
domain-wall valence quark masses as in the calculation of B^ presented in this work. We 
determine both and fx with ~ 3% accuracy, and their ratio with ~ 2% accuracy. Given 
the value of \Vud\ from superallowed /3-decay, our result for is consistent with experiment. 
Similarly, given the \Vus\ determination from semileptonic kaon decays using non-lattice 
theory, our result for fx is consistent with experiment. Our result for the ratio fK/f-w, which 
is independent of the CKM matrix elements, is consistent with other more precise lattice 



determinations 



24 



26|. Although our decay constant calculation does not check the Rome- 
Southampton nonperturbative renormalization (NPR) procedure, it does test the remaining 
ingredients in the calculation of Bk-, especially the chiral and continuum extrapolation using 
MA^PT. Therefore the successful calculation of the well-known quantities /tt and fx bolsters 
confidence in the calculation of the weak matrix element B^ presented in this work. 

This paper is organized as follows. In Sec. [Tll we describe the details of our numerical 
mixed-action lattice simulation; we present the actions and parameters used and describe 
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how the relevant 2-point and 3-point correlators are analyzed. Next, in Sec. Illlt we describe 
two independent calculations of the renormalization factor Zbj^ needed to match the lattice 
matrix element to the continuum. We compute using the nonperturbative Rome- 

Southampton approach; this is used in our preferred determination of B^. We also compute 
Zsji to 1-loop in mean-field improved lattice perturbation theory to provide a cross-check 
and aid in estimating the systematic error associated with the matching. In Sec. llVt we 
describe the extrapolation of Bk to the physical quark masses and the continuum using 
NLO MAxPT supplemented by higher-order analytic terms to allow an interpolation about 
the strange quark mass. Next, in Sec. IVj we present the systematic error budget for Bk, 
describing each individual uncertainty in a separate subsection for clarity. Finally, in Sec. IVIt 
we compare our results to previous unquenched lattice determinations and to the preferred 
values from the unitarity triangle analyses. We conclude by discussing the prospects for 
improvement in our mixed-action lattice calculation and for its phenomenological impact on 
the search for new physics. 



II. LATTICE CALCULATION 



In this section we describe the details of our numerical mixed-action lattice calculation. 
We first present the valence and sea quark lattice actions and input parameters (such as 
quark masses and lattice spacings) in Sec. Ill A[ Next, in Sec. Ill Bl we present the 2-point and 
3-point correlation functions needed to determine the unrenormalized lattice value of B^:- 
We give the valence quark propagator source wavefunctions and boundary conditions. We 
also describe the method used to extract Bk from a ratio of 3-point and 2-point functions, 
and show example correlated plateau fits with jackknife errors. 



A. Actions and input parameters 



We use the unquenched lattices generated by the MILC Collaboration for our numerical 
lattice calculation of Bk-, which include the effects of three dynamical flavors of asqtad- 
improved staggered fermions {2^]. Because the MILC configurations are available at several 
light quark masses and lattice spacings {g. I28I 



they allow us to have good control over 



the both the chiral extrapolation in the sea sector and the continuum extrapolation. We 
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TABLE I: Parameters of the MILC improved staggered gauge configurations and domain-wall 
valence quark propagators used to calculate the unrenormalized lattice value of Bk- Columns 
one and two list the approximate lattice spacings and lattice volumes (in lattice spacing units). 
Columns three and four show the nominal up/down {mi) and strange quark (m/i) masses in the 
sea, along with the corresponding pseudoscalar taste pion mass. Columns five and six list our 
partially quenched valence quark masses (mx), along with our lightest available domain- wall pion 
mass. Column seven shows the number of configurations analyzed on each ensemble. 



sea sector valence sector 



a(fm^ 


(-) 


ami/amh 


am-jt 


amx 


am-jt 


^conf. 


0.09 


40^ X 96 


0.0031/0.031 


0.10538(06) 


0.004, 0.0186, 0.046 


0.0999(12) 


150 


0.09 


28^ X 96 


0.0062/0.0186 


0.14619(14) 


0.0062, 0.0124, 0.0186, 0.046 


0.1212(17) 


160 


0.09 


28^ X 96 


0.0062/0.031 


0.14789(18) 


0.0062, 0.0124, 0.0186, 0.046 


0.1222(12) 


210 


0.09 


28^ X 96 


0.0124/0.031 


0.20635(18) 


0.0062, 0.0124, 0.0186, 0.046 


0.1216(11) 


198 


0.12 


24^ X 64 


0.005/0.05 


0.15971(20) 


0.007, 0.02, 0.03, 0.05, 0.065 


0.1718(11) 


216 


0.12 


20^ X 64 


0.007/0.05 


0.18891(20) 


0.01, 0.02, 0.03, 0.04, 0.05, 0.065 


0.1968(08) 


268 


0.12 


20^ X 64 


0.01/0.03 


0.22357(19) 


0.01, 0.02, 0.03, 0.05, 0.065 


0.1946(18) 


160 


0.12 


20^ X 64 


0.01/0.05 


0.22447(17) 


0.01, 0.02, 0.03, 0.05, 0.065 


0.1989(08) 


220 


0.12 


20^ X 64 


0.02/0.05 


0.31125(16) 


0.01, 0.03, 0.05, 0.065 


0.1949(13) 


117 



calculate on both the "coarse" (a ^ 0.12 fm) and "fine" (a ~ 0.09 fm) MILC ensembles, 
which have physical volumes ranging from approximately (2.5 - 3 fm)^. For each ensemble, 
the masses of the up and down sea quarks are degenerate; our lightest dynamical mass is 
approximately a tenth of the physical strange quark. For most of our ensembles, the mass of 
the dynamical strange quark is close to its physical value. At each lattice spacing, however, 
we have data on one ensemble with an unphysically light strange sea quark in order to better 
constrain the strange sea quark mass dependence and aid in the chiral extrapolation. The 
left-hand side of Table [T] shows the parameters of the MILC gauge configurations used to 
calculate Bk- 

We construct the 2-point and 3-point correlation functions needed to determine Bk us- 



ing domain-wall valence quark propagators 



29 



30| . The approximate chiral symmetry of 
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domain-wall fermions simplifies both the nonperturbative determination of the renormal- 
ization coefficient, ^_b^, and the extrapolation of Bk to physical quark masses and the 
continuum; these advantages will be discussed in greater detail in Sees. HII] and HVl respec- 
tively. We compute the domain-wall propagators using the Chroma software system for 
lattice QCD jsl. We use the same input parameters as the LHP Collaboration [sl; this 
allows us to check simple quantities such as the pion masses and the residual quark mass. 
We first smear the MILC lattices using the standard hypercubic blocking (HYP) parame- 
ters given in Ref. jsil in order to reduce the size of explicit chiral symmetry breaking and 
proximity to the Aoki phase 32|]. On both the coarse and fine ensembles we simulate with 
a domain-wall height of M^=l.l and a fifth dimension of length Ls=lQ. For each sea quark 
ensemble, we calculate at several valence quark masses; this allows us both to extrap- 
olate the numerical lattice data to the physical up/down quark mass and to interpolate to 
the physical strange quark mass. Our lightest valence quark mass is chosen to be as light 
as possible while keeping finite-volume effects under control. Specifically, we restrict the 
quantity m^L > 3.5 to keep 1-loop MA^PT finite volume effects for Bk below 1%. Thus 
the mass of our lightest domain-wall pion is ~ 280 MeV on the 2.5 fm ensembles, and ~ 240 
MeV on the 3.5 fm ensemble. The fifth column of Table|T]shows the bare domain-wall masses 
used to calculate Bk- 

In most mixed staggered sea, domain-wall valence lattice simulations, the bare domain- 
wall quark mass is tuned so that the mass of the domain-wall pion is equal to the mass of 
the lightest staggered pion in the sea sector jsj. Although this procedure does not eliminate 
unit arity- violating discretization effects in the mixed-action theory at nonzero lattice spac- 
ing, tuning the domain-wall pion to one of the staggered pion masses allows one to approach 
full QCD as the continuum limit is taken numerically, even for quantities for which mixed- 
action chiral perturbation theory expressions do not exist or are not applicable. Fortunately, 
for the case of Bk-, we can use MA^PT jol, to properly account for and remove these 
discretization errors in fits to quantities evaluated at multiple lattice spacings and valence 
and sea quark masses. Thus we do not make any attempt to tune the bare domain-wall 
quark masses in our lattice calculation. 

In order to convert dimensionful quantities calculated in our mixed-action lattice simu- 
lations into physical units, we need the value of the lattice spacing, a, which we determine 
by computing a known physical quantity that can be directly compared to experiment. Al- 
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though all of the coarse (or fine) MILC lattices have approximately the same lattice spacing, 
slight variations exist due to the choice of simulation parameters in the gauge action. We 
account for these differences by converting all of our data from lattice spacing units into 
ri units before perfor ming; any chiral fits. Because ri is related to the force between static 
quarks, rfF{ri] 



1.0 



33|, this method has the advantage t 



lat the ratio ri/a can be deter- 



mined precisely from a fit to the static quark potential [28|, |34]. The absolute scale, ri, can 
then be determined in several ways. In this work we use the scale ri = 0.3108(15) (lyg) fm 
to convert our simulation results into physical units. This value is obtained by combining 
the recent MILC determination of ri/^ with the experimentally measured value of /,r 24 1. 



We use an alternative determination of ri from the T splitting 



35 



ri = 0.318(7) fm, in 



order to estimate the systematic error due to the scale uncertainty. 



B. Three-point correlation functions 

Bk parametrizes the nonperturbative QCD contribution to CP-violating neutral kaon 
mixing. Kaon mixing occurs via electroweak box diagrams. Integrating out the heavy 
intermediate VT-bosons to isolate the hadronic contribution leads to the following = 2 
operator in the effective Hamiltonian: 

0^=' = [^7.(1 - 7M[S7,{1 - 7,)dl (1) 

where we omit the color indices for simplicity. In order to ensure that the value of is 
close to unity, Bk is defined ratio: 

R _ (K'\[s^,il-i,)d][s^,{l-^,)d]\K') 

= — =^ — — — -——^ K^) 



^{K |s7^(l -75)ci|0)(0|57^(l -75)c?|i^°) 
where the numerator is the desired AS = 2 matrix element, and the denominator is the 
same matrix element as in the numerator evaluated in the vacuum saturation approximation. 
Because the matrix element in the denominator is related to the kaon decay constant, Eq. ([2]) 
is often simplified as 

- « 2 f2 • l-^J 

3"''KJ K 

In this work we calculate Bk numerically from the following ratio of lattice correlation 
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functions: 



where T is the temporal extent of the lattice and we include the superscript "lat." to em- 
phasize that the quantity in Eq. (j4j) needs to be renormalized in order to recover Bk in a 
continuum regularization scheme. We fix the locations of the source and sink kaons in the 
numerator 3-point function at t^rc and t^^c + T, respectively, and vary the position of the 
four-quark operator, O^"^"^, over all time slices t in between. We use wall sources for our 
kaons throughout the calculation, but use local sinks for both the four-quark operator in 
the 3-point function and the axial-current operator in the 2-point functions: 

^p(*) = ^^i^^t)'y^d{y,t), (5) 

^Aif) = X]s(f,t)757^rf(f,t). (6) 

X 

The volume factor Li^ in the numerator of Eq. (jlj) accounts for the differing normalizations 
of the wall sources and point sources used in the determination of B^^\ 

For each domain-wall valence quark mass on a given MILC configuration, we compute 
two Coulomb gauge-fixed wall-source propagators starting from the same lattice timeslice, 
tgrc- one with periodic and another with antiperiodic boundary conditions in the temporal 
direction. The spatial boundary conditions are always periodic. The Coulomb gauge-fixed 
wall-source is used to reduce contamination from excited states. We then take symmetric 
and antisymmetric linear combinations in order to produce forward- and backward-moving 
propagators beginning at tsrc- We use these symmetrized propagators in the interpolating 
operators and (f) in order to effectively double the number of lattice timeslices. This 
ensures that finite-size effects due to pions circling the lattice in the temporal direction are 
negligible. Using the same time slice for the source of the forward- and backward-moving 
propagators also allows us to save a factor of two in computing time.^ 

In order to make the best use of our computing resources, we generate domain- wall quark 
propagators on every fourth recorded MILC gauge configuration (typically every 20th or 24th 
trajectory) in order to reduce autocorrelation errors. Our earliest runs have propagators with 



^ This method was suggested to us by N. Christ j37| . 
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tsrc = 0, which we chose for simphcity. In order to take advantage of the large temporal 
extent of the MILC lattices and further reduce autocorrelations, however, our later runs use 
a randomly chosen tgrc- Although the two data sets are expected to have somewhat different 
autocorrelation times, there is nothing a priori wrong with combining them in an ensemble 
average. 

Figured] shows a representative plateau fit on a coarse ensemble for 5^*/ (4L^) with a non- 
degenerate kaon made up of a light quark with mass around m^/G and a heavier quark with 
mass close to . Figure [2] shows a similar plateau fit on a fine ensemble where the heavier 
quark mass is again close to m^, and the light-quark mass is around nis/lO- The confidence 
levels of the fits are computed using the full correlation matrix in the minimization of in 
order to assess the quality of the plateaus. The statistical errors in the fit are determined by 
performing a separate fit to each single-elimination jackknife sample; the correlation matrix 
is remade for each jackknife fit. Excellent fits to a constant are found, and the confidence 
levels of the fits in Fig. [T] (CL = 0.71) and Fig. [2] (CL = 0.94) are typical of our numerical 
data. Although our plateau region appears to be quite long by inspection, a correlated 
fit requires a fit to a smaller span of the time extent so that the correlation matrix does 
not become too large to resolve with our current statistics (~ 150-270 configurations per 
ensemble). Thus, we limit our plateau fits to ~ 10-15 time slices. In practice, this is not 
much of a limitation, since we fold our data in the time direction. Typical statistical errors 
on the raw i?^* lattice data are at the sub-percent level, with 1-2% errors on the points with 
the lightest quark masses. 

Autocorrelation errors were studied on the two longest runs on the coarse {ami/anih = 
0.007/0.05) and fine {ami/am^ = 0.0031/0.031) ensembles. These errors are investigated by 
blocking the data before performing the single elimination jackknife estimate of the statistical 
error. However, there is also a correction to the statistical error coming from the fact that 
the correlation matrix is not known perfectly, but is determined approximately from the 
data set for a given fit. It has been shown in Ref. 38| that a jackknife fit with an estimate 
of the covariance matrix remade with each jackknife sample leads to a slight overestimate of 
the variance. This correction to the statistical error is at the ~ 5 — 10% level for our data set, 
and tends to cancel the expected (and difficult to resolve) slight increase in the statistical 
errors due to autocorrelations. Corrections to the statistical errors due to autocorrelations 
are at most a few percent. Given the rather small total correction to the statistical error 
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FIG. 1: Plateau fit to 5^7(4^3) on the coarse ami/arrih = 0.007/0.05 ensemble. The legend 
shows the non-degenerate pan of quark masses making up the kaon in the three-point correlation 
function. The correlated x^/dof and confidence level of the fit are given in the title. 
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FIG. 2: Same as Figured but on the fine ami/arrih = 0.0031/0.031 ensemble. 
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due to the combination of these effects, we do not adjust the errors "by hand" as input to 
later (chiral and continuum extrapolation) fits. This correction to the statistical error due 
is also a small fraction of our final total error, and can be neglected. 



III. RENORMALIZATION OF THE AS = 2 OPERATOR 

In this section we describe the calculation of the renormalization factor, Zb^, which 
is needed to match the lattice matrix element to the continuum. We present the result 
renormalized in the MS scheme at 2 GeV. We determine using two independent meth- 
ods: lattice perturbation theory and the nonperturbative Rome-Southampton approach. 
Although we use the nonperturbatively determined Zbj^ to calculate our central value for 
Bx, the lattice perturbation theory calculation provides a valuable crosscheck on the non- 
perturbative renormalization and an indication of the size of the systematic uncertainty on 
the renormalization factor. 



A. Lattice perturbation theory calculation of Zbj^ 

In this subsection, we use lattice perturbation theory to match our lattice calculation of 
Bk to the MS scheme using naive dimensional regularization (NDR). Although naive lattice 
perturbation theory appears to converge slowly, two main causes of this have been identified 
in Ref. 39|. The first is that the bare gauge coupling is a poor expansion parameter, and the 
second is that tadpole diagrams tend to have large coefficients. If a renormalized coupling is 
used and one restricts oneself to quantities for which tadpole diagrams largely cancel, then 
lattice perturbation theory appears to converge as well as continuum perturbation theory. 
The difficulties with large tadpole corrections are present even in chiral fermion formulations, 
where they are just as serious as in other formulations. We address this issue here. 

In the case of domain-wall quarks, the domain-wall height M5 is additively renormalized, 
and large tadpole corrections can appear. It has been shown that mean-field improvement 
is then necessary to get correct results from lattice perturbation theory for domain-wall 



quarks |40l442l] . Our calculations do not suffer from large tadpole corrections because we 
use HYP-smeared domain-wall quarks in our simulations. This HYP-smearing smoothes the 
gauge fields, and has the effect of dramatically reducing the tadpole contributions in lattice 
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FIG. 3: Vertex diagram for the correction to bilinear operators in lattice perturbation theory. 



perturbation theory for our simulation parameters. Thus, the difference between naive and 
mean-field improved lattice perturbation theory in our renormalization of Bk is small. Even 
so, we adopt the correct mean-field improvement in all results presented here. 

The renormalization factor matching the lattice calculation of Bk to the MS scheme can 



be written as 43 1 



(1-wg) ^ZjZ+{fia) ^ Z+{fia) 

where Z+ is the renormalization factor for the operator O^^^^, Za renormalizes the axial 
current, Wq = 1 — M5, and Z^ is the quantum correction to the normalization factor 1 —w^. 
It is useful to define in this way, since the tadpole and self-energy corrections cancel. 
The renormalization factor contains the running of the operator from the lattice scale to 
the continuum scale /i. The relevant Feynman diagram for our particular lattice calculation 
is shown in Fig. [3] and the necessary Feynman rules are given in Appendix \M In the MS 
scheme with NDR, we have {43 1 



-41n(/ia) - H + 2 Invr^ + '^{167r^){Is 



(8) 



with Isy defined in Eq. (IA37p . 

Given the cancellation of tadpoles in Eq. ([7]), the only effect of mean- field improvement 
in the one-loop renormalization factor Zbj^ is to shift the domain-wall height M5 — = 
M5 — 4(1 — uq), where uq is the fourth-root of the plaquette. This shift would be appreciable 
if not for the HYP smearing of the domain wall quarks, since Uq ~ 0.87 on the MILC coarse 
and fine lattices. However, for our calculation it is appropriate to take the HYP-smeared 

13 



plaquette in the mean-field improvement factor, and tliis leads to y^^''^"'^'^^'^ = 0.984 and 
^MF,finc _ g ggY g^-^^ ^ significantly smaller shift in . The final result for decreases 
by only about one percent (at both coarse and fine lattice spacings) after adopting the 
mean- field improvement. 

We adopt the Brodsky, Lepage, and Mackenzie (BLM) scheme for setting the scale in the 



coupling Oiy[g{q*) [4^. The BLM prescription for obtaining q* is 



where f{k) is the one-loop integrand, and the numerator is the first log moment. Note that 
throughout this section and in Appendix |Xl all momentum integrals run over — tt < < it, 
with kn in lattice units. This prescription for computing q* is well-defined for finite lattice 
integrals. In the case of B^, however, where there is an anomalous dimension, the BLM 
pres crip tion needs to be modified. We follow the prescription introduced by C. Bernard et 



al. 



45| . and discussed in detail by DeGrand 



46| . A generic integral evaluated in the MS 



scheme will take the following form: 



, ^2a; u 1 



A^i-7£ + log(4ir)|+^log^ + 4 + S, (10) 



where 2a; = 4 — 2e is the dimension of the integral and where the term in curly brackets is 
discarded to give a finite integral, J^. The log moment of the divergent part of the one- 
loop expression must be handled with care. The log moments of the finite lattice integrals 
ly and Is are straightforward to evaluate using Eq. ([9]), and are denoted ly and Ig in 
Table Ull However, we need the log moment corresponding to the entire term in brackets 
in Eq. (|8]), including the anomalous dimension. The log moment corresponding to the 
anomalous dimension [the first three terms in brackets in Eq. (|H])] can be defined as the log 



moment of the following finite integral 46 1 



I^ = AJ, + Bh, (11) 



We actually compute av{q*): the strong coupling constant in the V scheme, and exploit the fact that 
'^MS ~ '^^ order we are working. The scales used to determine the coupling in each scheme are 

related by ln[(agy)^] = ln[(agijg)^] + 5/3 in the first order BLM prescription where onfj^ the first log 



moment is required. The V scheme is defined with respect to the heavy-quark potential 39l. |44 1 
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TABLE II: Computed values of Z^j^ in the BLM prescription. The first column labels the approx- 
imate lattice spacing in fm. The second column is the numerical evaluation of the integral ly, and 
the third is that of the integral !$■ The fourth and fifth columns are the first moments of ly and 
Is, respectively. The sixth column is o-Qblm^ ^"^^ seventh column is Z^^'^^^(2GeV). Errors 
from numerical approximation of the integrals are no more than one digit in the last displayed 
decimal. 



a (fm) 


Iv 


Is 


T* 


T* 
^S 




^MS,NDR^2GeV) 


0.12 


0.0158 


-0.0161 


0.0336 


-0.0150 


1.56 


0.909 


0.09 


0.0158 


-0.0155 


0.0336 


-0.0154 


1.42 


0.955 



where the F stands for finite, and 



Jl 

J2 



167r^ 



167r' 



(2^ 



1 



(fc2)2 
1 



(P + /i2)2_ 
1 



P(A;2 + /i2) (F + /i2)2 



(12) 
(13) 



where /i is the MS scale (in lattice units), and 9 is the Heaviside step function. The values of 
A and B for Zbj^ in the MS, NDR scheme are —2 and —5/3, respectively. The log moment 
of the one-loop expression for Zbj^ can then be used to compute q* using 

(^&)* + i(16vr^)(^5-^y) 



*\2 



ln(ag*) 



—4 In(yua) 



11 

3 



21n7r2 + |(167r2)(/5-/v 



(14) 



where (-^^)* signifies that the first log moment is taken in the momentum integrals appearing 
in Eq. ( fTTl) . The computed values for the integrals ly and J5, as well as their first log 
moments, are given in Table HTl The resulting g*'s and the final values for z^'^^^^ are also 
given. All integrals were evaluated numerically using the Mathematica package and 



results in Ref. 



43| using the same action but without HYP-smearing were reproduced. 



B. Nonperturbative determination of Zbi, 



1. Rome- Southampton method 



We compute the renormalization coefficient for Bk nonperturbatively in the RI/MOM 



scheme devised by the Rome-Southampton group 



10|. In this scheme, the simple renormal- 
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ization condition is that the renormahzed n-point functions in Landau gauge are equal to 
their tree-level values. Because the RI/MOM scheme is regularization-invariant, it is useful 
for both perturbative or non-perturbative techniques. Thus it is well-suited to lattice QCD 
simulations. Once Zbj^ has been determined nonperturbatively in the RI/MOM scheme, it 
can easily be converted to the MS scheme and run to the scale /i = 2 GeV using continuum 
perturbation theory. 

The Rome- Southampton nonperturbative renormalization technique has already been 
successfully applied to lattice QCD calculations of with domain wall valence and sea 



quarks by the RBC and UKQCD Collaborations in Refs. (48l . l49l |. We can determine the 
renormalization factor Zbj^ in the same simple manner for our mixed-action lattice QCD 
simulations because the properties of the mixing coefficients are largely determined by the 
symmetries of the valence sector. In particular, errors of 0(a) are suppressed by ~ e~"'^^ 
Furthermore, mixings between the desired four-fermion operator and other operators of 
incorrect chirality are suppressed due to the approximate chiral symmetry, as we show in 
Appendix O 

Our primary goal is to determine the renormalization coefficient of the four-quark oper- 
ator given in Eq. ([1]): 

where we now show explicitly the chirality of the operator. Because chiral symmetry is 
slightly broken in our simulations, however, this operator can mix with other AS* = 2 
operators of different chiralities: 

Ovv-AA = [s7m(1 - 15)(^ [s7m(1 + l5)d], (15) 
Oss±pp = m - 75)dMl T 75)^^], (16) 
Ott = [sa^uil - l5)d] [sa^uil - 75)^]- (17) 

Thus the renormahzed B^ operator in principle receives contributions from all of the above 
operators, and is given in terms of bare lattice operators by 

C'x" = X] ^VV+AA,i Oi , (18) 

i 

where i E {VV + AA, VV - AA, SS + PP, SS - PP, TT}. The theoretical arguments of 
Ref. |49| suggest that the wrong chirality mixing coefficients should be of 0[{amj-es)'^], and 
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our data is consistent with this expectation within statistical errors. Thus, despite the fact 
that chiral perturbation theory predicts the corresponding B parameters of these operators 
to diverge in the chiral limit 48, SqI, their contributions to are negligible. 



In order to calculate via the Rome-Southampton approach 10|], we first compute 
the 5x5 matrix 

M,, = p,[rl^**], (19) 

where F-*^** is the amputated four-point Green's function in momentum space, i,j G {VV + 
AA, VV — AA, SS + PP, SS — PP, TT}, and the projector Pj selects out the component wit 



chirality j. The details of this procedure are described thoroughly in Appendix B of Ref. [48] . 
We also compute the tree-level value of this matrix by setting all of the momentum-space 
propagators in the amputated Green's functions equal to the identity: 

^tree _p.^ptree|_ (20) 

We then impose the RI/MOM renormalization condition, 

^M,. = Mr , (21) 
where Zg is the quark wavefunction renormalization factor, in order to obtain the quantity 

^ = Mr M-/ . (22) 



The renormalization coefficients for the various four-fermion operators are then given by 

(23) 



7 7 / 7 ^ ^ 



zl zi 

where is the renormalization factor for the axial current. For example, the dominant 
contribution to the Bk lattice operator renormalization comes from the diagonal mixing 
coefficient: 

Z,,„ = = . (24) 



In order to determine the four-fermion operator mixing coefficients using Eq. fl23|l . we 
also need the ratio Zq/Z^- Fortunately, the renormalization factors for the quark bilinears 
can also be calculated in a simple manner using the Rome-Southampton method. The 
renormalization coefficients relate the bare and renormalized quark bilinear operators in the 
following manner: 

[uFcijren = ZY\jlTd]Q . (25) 
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In order to determine Zr, we first compute the bare Green's functions between off-shell 
quarks in momentum space. We then amputate the Green's function and separately project 
out the components with each chirality to obtain the bare vertex amplitudes: 

1 



A5 
Ap 

Ay 

Aa 
At 



12 
1 

12 
1 

48 



Tr[Gr''l] , 
Tr[G7P75] 



Tr 



It, 

48 



It, 

72 



Eg: 
Eg: 



amp 

7m 



7a. 



7m 75 



amp 
M-v '^M 



(26) 
(27) 

(28) 
(29) 
(30) 



Finally, we impose the RI/MOM renormalization condition 



A,. 



-A,, 



(31) 



(32) 



The renormalization coefficients for the quark bilinears are then given by 

Zg Ai' 

In the RI/MOM prescription, the four-fermion operator renormalization coefficients are 
given as functions of the momenta in the amputated Green's functions used to determine 
Zij/Zg, which are chosen to be identical for all four quarks in our computation of Zbj^- 
We therefore need to extract at a sufficiently high momentum that hadronic effects 
are negligible and the momentum-dependence can be described by perturbation theory. We 
cannot use too high a momentum, however, or lattice discretization errors will become large. 
Thus use of the Rome-Southampton technique requires the existence of a momentum window 
in which both hadronic effects and discretization errors can be neglected: 



Aqcd < J5 < a" 



(33) 



In practice, however, we need to work in the region (ap)"^ > 1 in order to avoid large violations 
of chiral symmetry which we observe at low momenta. Fortunately, discretization effects in 
the region of interest, p ^ 2 GeV, are generally rather small and can be taken into account 
by a simple linear fit in (ap)^, as discussed in Ref. 51|- This is the approach that we take 
in the calculation of in Sec. IIIIB 31 
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2. Chiral symmetry breaking and — Ay 



Although the calculation of Zb^ only requires the renormalization factor for the axial 
current, Z^, the vector and axial- vector current renormalization factors should be equal in 
the chiral limit for sufficiently large momenta due to chiral symmetry: 

Za = Zv , (34) 

or equivalently, 

Aa = Av . (35) 

Thus we can take the average of these two quantities in order to reduce the statistical error 
in Zq/ZA using the relationship 

Aa = ^{Aa + Av) ■ (36) 

In practice, however, ^ Ay in the chiral limit for any value of the momentum in our 
nonperturbative determination. Figure H] shows the extrapolation of the quantity 2{Aa — 
Ay)/(A^ + Ay) on the coarse lattice to the chiral limit at p ^ 2 GeV using a function that 
is linear in both the valence and sea quark masses. This quantity provides an indication 
of the amount of chiral symmetry breaking in the computation. At nonzero quark mass, 
2(Aa— Ay) / (A^+Ay) can be as large as ~ 1% in the momentum region (ap)"^ > 1 that we are 
using to extract Zbj^- The difference between A^ and Ay decreases towards the chiral limit, 
as is expected, but never becomes consistent with zero. Figure [5] shows 2(A^— Ay)/(A^+Ay) 
versus (ap)^ on the ami/am^ = 0.007/0.05 coarse ensemble for the five available valence 
quark masses and in the chiral limit. Again, it decreases in magnitude as expected at larger 
momenta, but is never zero. We observe similar behavior on the fine lattice. This persistent 



difference between A^ and A 
UKQCD collaborations 
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las also been observed and studied in detail by the RBC and 



51 



52| . and can be attributed to several sources. 
The first source is explicit chiral symmetry breaking due to the nonzero quark masses 
used in simulations 49|. Use of the operator product expansion shows that A^ and Ay can 
receive contributions proportional to 

mg{qq) 



p2 ' p4 



(37) 



at lowest order in l/p"^ [51|. Because these operators are proportional to rrig, they explicitly 
break chiral symmetry and need not contribute equally to A^ and Ay. The contribution of 
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FIG. 4: Chiral extrapolation of 2(A^ — Ay)/(A^-|-Ay) on the coarse lattice at (ap)'^ = 1.468 using 
a linear function in nix and nii. Although only the data points with filled symbols were used in the 
fit, the fit line does a reasonable job of describing the heavier data points that were not included. 
The cyan error band shows the extrapolation/interpolation for points where the domain- wall pion 
mass is tuned to equal the lightest (taste pseudoscalar) staggered pion mass. 



the operators in Eq. (137|) can be seen clearly in the data. As Fig. [5] shows, the difference 
between and Ay increases rapidly as the momentum approaches zero and decreases slowly 
as the momentum becomes larger than ^ 2 GeV. The contributions from the operators in 
Eq. dnZD can be removed by extrapolating to the chiral limit at fixed momentum. Figs. H] 
and [5] show that, although this procedure does indeed reduce the splitting between A^ and 
Ay, it does not eliminate the difference. Thus there must be additional sources of chiral 
symmetry breaking, which we now discuss. 

The next source is chiral symmetry breaking due to the use of a finite Lg. Theoreti- 
cal arguments suggest, however, that this would lead to errors that are much smaller, of 
(9((ami.cs)^) ~ 10^^ in our numerical simulations 
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52l |. This would produce a negligible 



difference between A^i and Ay, and cannot account for the size of the difference that we 
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FIG. 5: The quantity 2{Aa — Ay)/{Aa + Ay) versus {ap)^ on the ami/amh = 0.007/0.05 coarse 
ensemble for several valence quark masses and in the chiral limit arrix = ami = 0. 



observe in the data. 

A more significant source of chiral symmetry breaking that does not vanish in the chiral 
limit is the choice of kinematics used to compute both Zij/Z^ and Aj. As in the standard 
RI/MOM prescription, we are using "exceptional momenta" configurations in which there 
is no momentum transferred between the initial and final states, or more precisely 

Pi=Pf=P , (38) 



where Pij are the momenta of the initial and final states. It was shown in Ref. 49|] that 
this, unfortunately, leads to contributions to A^ — Ay of the form 

(qqqq) 



p2 



(39) 



Because this operator is not proportional to the quark mass, it does not vanish in the chiral 
limit at fixed momentum. This contribution can be removed by performing the nonpertur- 
bative renormalization calculation at non-exceptional kinematics, in which the sum of any 
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FIG. 6: The quantity 2(A/i— Ay) /(A^+Ay) versus (ap)"^ computed using non-exceptional momenta 
on the ami/anih = 0.007/0.05 coarse ensemble for several valence quark masses and in the chiral 
limit arux = ami = 0. 



subset of external momenta is nonzero. In this case we have 

P^i =P} = {Pi-PfY =P\ (40) 

but Pi ^ Pf- In order to check that this is indeed the source of the difference between 
Aa and Ay in our data, we have also computed 2{Aa — Av)/{Aa + Ay) at non-exceptional 
kinematics; the results are shown in Fig. [6l Although the statistical errors are not as small, 
the results in the chiral limit are consistent with zero for sufficiently large values of (ap)"^. 

For the calculation of Z^^, in this work, we use exceptional kinematics, despite the result- 
ing chiral symmetry breaking. This is because the continuum perturbation theory needed to 
convert the result from the RI/MOM scheme to the MS scheme has not yet been calculated 
for non-exceptional kinematics.^ We therefore include the difference between A^ and Ay as 

^ The expressions needed to convert the quark bilinears from the RI/MOM scheme to the MS scheme have 
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a source of systematic uncertainty in Zbj^- 



3. Nonperturbative renormalization factor calculation 

We now present the nonperturbative determination of Zbi^, which we compute from the 
quantity 

Zvv+AA,VV+AA I + ^^^^ 



Zl \ 2 

using the method of Rome-Southampton. Table UTTl shows the parameters used in generating 
the NPR lattice data set. We have several valence and sea quark mass combinations on both 
the coarse and fine lattices in order to allow an extrapolation of Zb^^ to the chiral limit. For 
those ensembles that are listed as "blocked" in the table, we computed Zb^ on every sixth 
trajectory and blocked the data in groups of four in order to reduce autocorrelations. On 
those ensembles for which the data was not blocked, we computed Zbj^ only on every 24*^^ 
trajectory. 

Because we must extrapolate Zbj^ to the chiral limit in both the valence and sea quark 
sectors, on the coarse lattice we have generated data on three ensembles at the nominal 
strange quark mass {artih = 0.05) and on one ensemble with a lighter than physical strange 
sea quark mass (am^ = 0.03). At our current level of statistics the results for Zbj^ on the 
ami/am^ = 0.01/0.05 and ami/arrih = 0.01/0.03 coarse ensembles are consistent. Because 
we do not observe any dependence on the strange sea quark mass in our data, we fit our data 
assuming only a dependence on the light sea quark mass to determine the central value of 
Zbk- We use an alternative fit that includes strange sea quark mass dependence to estimate 
the systematic error associated with extrapolating the strange sea quark mass to the chiral 
limit. 

We first extrapolate Zb^ using a polynomial function in the valence and sea quark masses: 

/nsea,n.al (Xval, Xsea; P^) = Zf^^°^{p'^) + ^ Ci,sea(P^)Xsea + <^i.val(p^)Xval . (42) 

1=1 1=1 



recently been computed to one-loop order for non-exceptional kinematics by Sturm et al. in Ref. 53 1 
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TABLE III: Lattice QCD data used in the nonperturbative renormalization of B^- For those 
configurations that were blocked, Zbj^ was computed on every 6^^ configuration and blocked in 
groups of 4. For those configurations that were not blocked, Zbj^ was computed on every 24*^^ 
configuration. 



a(fm) (^) X ^ ami/arrih anix -^conf. blocked? 



0.09 


28^ X 96 


0.0062/0.031 


0.0119,0.0171,0.0287,0.04 


387 


no 


0.09 


28^ X 96 


0.0093/0.031 


0.0287 


251 


no 


0.09 


28^ X 96 


0.0124/0.031 


0.0287 


381 


no 


0.12 


20^ X 64 


0.007/0.05 


0.01,0.02,0.033,0.038,0.056 


836 


yes 


0.12 


20^ X 64 


0.01/0.05 


0.01,0.02,0.033,0.038,0.056 


540 


yes 


0.12 


20^ X 64 


0.02/0.05 


0.01,0.02,0.033,0.038,0.056 


484 


yes 


0.12 


20^ X 64 


0.01/0.03 


0.01,0.02,0.03 


81 


no 



where 

Xval = ^^yy[2(m^+mres)] , (44) 

are dimensionless ratios < 1, and the parameters /istag and 

/Udw are obtained from tree-level 
MAxPT fits to the pseudoscalar meson masses. In order to determine the preferred fit 
ansatz, we independently increase nsea and riyai until the correlated confidence level of the fit 
no longer increases significantly. We find that this occurs when n^ca = 1 aud n^ai = 2. We 
use the fits with additional terms to estimate the systematic uncertainty due to the choice 
of fit function. 

Figures [7] and [8] show the chiral extrapolation of Zbj^ on the coarse and fine lattices 
at (2 GeV)^ using a fit function linear in the light sea quark mass and quadratic in 

the valence quark mass. Although this extrapolation is nominally at "fixed (ap)^", this is 
not quite true. This is because, although all of the coarse (or fine) MILC ensembles are 
generated with approximately the same lattice spacing, there are slight fluctuations in the 
lattice spacing from ensemble to ensemble. Thus data on different ensembles with the same 
value of {apY do not correspond to precisely the same physical momentum. We convert our 
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FIG. 7: Chiral extrapolation of -^B// coarse lattice at (ap)^ = 1.468. Note that the 

fit lines for the ami/am^ = 0.01/0.05 and ami/am^ = 0.01/0.03 ensembles are identical because 
we have not included any strange sea quark mass dependence in the fit function. The cyan band 
shows the extrapolation along the trajectory m™^ = m^"'. 



data into ri units using the value of ri/a determined on each ensemble before performing 
the chiral extrapolation, so that everything is in the same units. Fortunately, the variation 
in ri leads to only a slight variation in the momentum- squared, of ~ 0.1%. Because this is 
even smaller than the statistical errors in our data points, the resulting systematic error can 
be safely neglected. 

Next we attempt to remove discretization errors in due to the fact that we are 
extracting Zbj^ at momenta that are of 0{a~^). Following the same procedure as in Ref. 49| . 
we use the continuum 1-loop perturbation theory expressions in Eqs. (IB5p - (lB8p to convert 
^bI^^*^^ to Z^^. This is shown in Fig. P (Fig. [TU]) for the coarse (fine) lattice. Although the 
quantity Zf^^ should be scale-invariant, we observe that Z^^ in fact has an approximately 
linear dependence upon (ap)"^ in the region of interest. We believe that this scale-dependence 
is primarily from lattice artifacts that can be removed by performing a linear extrapolation 
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C.L. = 0.35 




shows the extrapolation along the trajectory = m: 



in (ap)^ . We therefore fit the data to the form 

A + B{apf (45) 

and interpret the intercept A as the true Z^^. We extrapolate to its true value as 
shown in Figs. [3 and [10] using the momentum range 2 GeV ^ p ^ 2.5 GeV. This choice 
satisfies the criterion p ^ Aqcd needed to avoid hadronic effects. Specifically, for the coarse 
data, we fit within 1.5 < (ap)"^ < 2.3 and for the fine data we fit within 0.8 < (ap)^ < 1.2. 
Variations of these fit regions do not alter the final results significantly. We obtain 

^Sl^oarse ^ 1,2822(29) , 

^sjMe ^ 1.3033(93) , 

where the errors are statistical only. We note, however, that some of the scale-dependence 
in Z^^ may in fact be due to the lack of higher-order terms in the matching factor. We 
therefore account for this and other errors due to the truncation of perturbation theory in 
the systematic error budget for Zb^^- 
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(ap)' 



FIG. 9: Linear-in-(ap)^ extrapolation of on the coarse lattice to remove lattice artifacts using 
the fit range 1.5 < {ap)'^ < 2.3. The data points used in the fit are denoted by the filled circles and 
the true value of Z^^ obtained at {ap)"^ = is denoted by the star. 

Finally, we convert Z^^^ to Z^^ and run it to 2 GeV again using Eqs. f lBSp - flBSp . We 
obtain 

ZU'™(2 GeV) = 0.9339(21) , 
^^,fino^2 GeV) = 0.9493(68) , 

where the errors are statistical only. We estimate the systematic uncertainties in 
Z^{2 GeV) later in Sec. |VDl 

IV. DETERMINATION OF Bk 

In this section we describe the extrapolation of Bk to physical quark masses and the 
continuum. In Sec. IIV At we present the expression for Bk at next-to-leading order (NLO) 
in MA^PT and describe those features that are most relevant for the chiral-continuum 
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FIG. 10: Linear-in-(ap)^ extrapolation of on the fine lattice to remove lattice artifacts using 
the fit range 0.8 < (ap)'^ < 1.2. The data points used in the fit are denoted by the filled circles and 
the true value of obtained at {ap)^ = is denoted by the star. 



extrapolation. We then discuss the details of the chiral-continuum extrapolation procedure 
in Sec. HVBl 



A. Bk at NLO in MAxPT 

We first review the tree-level mass relations for light pseudoscalar mesons in MA^PT since 
they are useful in understan ding the leading-order lattice-spacing contributions to mixed- 



action numerical simulations 22|. In a mixed-action theory one can have mesons composed 
of two sea quarks, two valence quarks, or one of each. At tree-level in MAxPT, discretization 
effects lead to different additive shifts to the masses of the three types of mesons. These 
mass-shifts are the only new parameters as compared to the continuum at this order, and 
their values have all been determined for our choice of mixed-action simulation parameters. 
The tree-level mass-shifts on both the coarse and fine MILC lattices are given in Table IIV[ 
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TABLE IV: Tree-level mass-shifts on the coarse and fine MILC lattices for our choice of mixed- 
action simulation parameters. The taste-singlet mass-splitting, A/, is independent of the valence 
sector [7|, while the residual quark mass, nij-es, and mixed meson mass-splitting, A^ix, both depend 
upon the choice of HYP-smearing • Errors shown are statistical only. 



a(fm) 


sea sector 


valence sector 

rimres 


mixed sector 
rja^ Ajnix 


0.12 
0.09 


0.537(15) 
0.206(17) 


0.0044(1) 
0.0016(2) 


0.207(16) 
0.095(20) 



In the sea sector of the mixed action theory, each flavor of staggered quark comes in 
four species, or "tastes" ; consequently, each flavor of staggered pseudoscalar meson comes in 
sixteen tastes. In the continuum, these tastes are identical and are related by an SU (4) sym- 
metry At nonzero lattice spacing, however, discretization effects split the degeneracies 



among the sixteen pseudoscalar meson tastes 541 ]: 



mss',t = /^stag("^s + my) + a At, (46) 



where s and s' are the staggered quark flavors, /istag is a regularization-dependent continuum 
low-energy constant, and Aj is the mass-splitting of a pion with taste t. At leading-order in 
staggered xPT (S^PT), a residual 5*0(4) taste symmetry ensures that the mass-splittings 



are identical for pions in the same S'0(4)-irrep: P, V, A, T, / [5J]. An exact U {1)a symmetry 
protects the taste pseudoscalar meson from receiving a mass-shift to all-orders in S^PT, 
implying that Ap = 0. 

At NLO in the mixed-action theory, the only nonzero staggered mass-splitting that is 



relevant is that of the taste-singlet, Aj 22|]. This is because the domain- wall valence quarks 
do not carry the taste quantum number; therefore mixed valence-sea four-fermion operators 
must contain two domain-wall quarks and two taste-singlet staggered quarks in order to be 
overall taste- invariant. The mass-splitting Aj has been calculated by the MILC Collabora- 
tion on both the coarse and fine MILC lattices [3], and is given in Table IIVI Because the 
parameter A/ is known, we reduce the number of free parameters in the chiral and continuum 
extrapolation of Bx by fixing Aj to the values in Table IIVI The mass-splitting Aj turns out 
to be the largest of the taste-splittings, and comparable to the taste-Goldstone pion mass on 
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the coarse MILC lattices, so the taste-singlet sea-sea mesons are quite heavy on the coarse 
lattices. Because the mass-splittings arise from discretization effects, however, they become 
smaller as the lattice spacing decreases. Specifically, the staggered taste-splittings scale as 
0{a'^a'^) since the Asqtad staggered action is (9 (a^)- improved. Thus a^A/ is already more 
than a factor of two smaller on the fine MILC lattices than on the coarse. 



In the valence sector of the mixed-action theory, domain-wall quar 



<;s receive an additive 



contribution to their mass from explicit chiral symmetry breaking 
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30|: 



m^v' = l^dwirriy + m^/ + 2mi.es), (47) 

where v and v' are the domain-wall quark flavors and mros is the residual quark mass. The 
size of rrires parametrizes the amount of chiral symmetry breaking in the valence sector, and 
is controlled by the length of the fifth dimension. We have determined the value of m^es in 



our mixed-action simulations in a previous publication (Ref. [17|) and present the results in 
Table IIVI On the coarse MILC lattices, we find that the value of mres in the chiral limit 
is about 3/4 the physical light quark mass; rrij-cs is approximately a factor of three smaller 
on the fine MILC lattices, i.e. 1/4 the physical light quark mass. The small values of the 
residual quark mass indicate that chiral symmetry breaking is under control in our mixed 
action lattice simulations. 

In order to reduce the number of fit parameters in our chiral and continuum extrapolation 
of Bk, we fix the value of mres in our chiral fits. We do not, however, use the values of mres 
given in Table llVt which are found by taking the chiral limit (m/ = ruh = rrix = 0) in both 
the valence and sea sectors. Instead, for each lattice data point, we fix mres to the value 
determined at that particular combination of valence and sea-quark masses. This effectively 
includes higher order corrections to mres and improves the confidence levels of our chiral fits 
to and m^. 

Because the mixed-action lattice theory has new four-fermion operators, the chiral effec- 
tive theory has new low-energy constants due to discretization effects. It turns out, however, 



that the mixed-action chiral Lagrangian has only one new constant at lowest order [22|. This 
coefficient combines with coefficients coming from the taste-symmetry breaking operators in 
the staggered sector SsJ to produce an 0{a'^) shift to the mixed valence-sea meson mass- 
squared: 

ml^ = fJ^dwim^ + m^es) + /istag"^s + a^Amix, (48) 
30 



ated the value of Amix for the parameters of our mixed-action lattice simulations in 
l?! . and present the results in Table HVl We find that the size of A^ix is less than half 



where v is the domain-wall quark flavor, s is the staggered quark flavor, and Amix is the 
effective mixed valence-sea meson mass-splitting obtained in lattice simulations.^ We have 
calcu 
Ref. 

of the taste-singlet staggered mass-splitting. A/, on both the coarse and fine MILC lattices. 

We do not need to fix the value of Amix during the chiral and continuum extrapolation 
of Bk because it turns out that the parameter Amix does not enter the expression for Bk 
in MAxPT at NLO, Eq. (149|) Although the mass-splitting enters the mixed-action 
expression for fx, it cancels exactly at NLO between the numerator and denominator in the 
ratio of matrix elements that defines B^, Eq. ([2]). 

Finally, we note that, for the purpose of our chiral and continuum extrapolation of B^, 
it is useful to express the tree-level meson masses in terms of the bare lattice quark masses 
given in Table HI not in terms of the renormalized quark masses. Because the valence and 
sea quarks are renormalized according to different schemes, we absorb the scheme-dependent 
quark-mass renormalization factors into separate coefficients of proportionality, /Xdw and 
/istag, in the tree-level mass relations, Eqs. (l46l)-(H8ll. 

The NLO xPT expression for B^ in a mixed-action domain-wall valence, staggered sea 



theory with 2+1 flavors of dynamical sea quarks is {q] 

p N PQ,2+1 . 

^] - 1 + . [/ + /2+11 

bJ ~ ^ 167r2/.>^, [^conn + l,,,,\ 

2 ^2 \2 



J X 



(49) 

'It]" 

I xy L xy 

where t^\(y) mass-squared of a meson composed of two x{y) valence quarks and 

fn'ipl^Hp) mass-squared of a taste-pseudoscalar meson composed of two l{h) sea quarks. 

The 1-loop chiral logarithms are separated into contributions from quark-level connected 
and disconnected diagrams. The parameter Bq is the tree-level value of Bk obtained in the 
continuum and SU{3) chiral limits. The four analytic terms, Ci - C4, are the only additional 
free parameters in the expression for Bk at NLO. Although the analytic term proportional 



^ The explicit expression for the effective mixed meson mass sphtting Amix, given as a Unear combination 
of the staggered sea taste-splittings and the new splitting unique to the mixed action theory, is derived 



in Ref. 55 
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to 0? is not present in the continuum, it is present for chiral lattice fermions. Thus Bk in 
the mixed-action theory, Eq. f H9|) . has no more undetermined coefficients than in the purely 
domain- wall case. 

The connected contribution to Bk is 

I conn = ^m\~l{m%) - i{m\){m\ + m^^) - £(m^)(m^ + m^^). (50) 

The chiral logarithms, i and are defined as 

01 2^ 2 A , .FV( .FVi T^ 4 ^ Ki{\r\mL) 

^(m ) = m ( In— + 5i (mL) 1 , 5^ ("^L) = — , (51) 

= -(ln^ + l)+53^^(mL), 53^^(mL) = 2 Ko(|r>L) , (52) 

where the difference between the finite and infinite volume result is given by 5f^(mL), 
and Kq and Ki are modified Bessel functions of imaginary argument. The disconnected 
contribution to Bk is 



I^^sl - 3« 5^2^ 5^2^ 



\^i{m]) {ml^ + m'^Rf'\{M%]■{^^f])^, (53) 



where 



To] 

WxY,i} = {mx,mY,mnj}, (55) 
{/if} = {mL^^niHj}. (56) 

When the up and down sea quark masses are degenerate, the fiavor-neutral, taste-singlet 
mass eigenstates are 

< = ^ + (57) 
and the disconnected contribution to Bk simplifies: 

Iltsl = 1{Ix + Iy + Ir,) , (58) 
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Iy = ix{x ^ r), 



(60) 



£(m: 



2 ^2 



(61) 



All of the sea quark dependence in the chiral logarithms appears in the disconnected terms, 
the sum of which vanishes for degenerate valence quark masses. The contribution vanishes 
identically when mx = my- In the limit that mx — >■ "ly, Ix — )■ —Iy- Thus the sum 
Ix + Iy + Iv^ 0. 



B. Chiral and continuum extrapolation of Bk 

We use the SU{3) MA^PT formula of Eq. ( H9ll in the extrapolation of our numerical 
lattice data to the continuum and to physical quark masses. The choice of SU{3) xPT 
is appropriate given the parameters of our numerical simulations because our light pion 
masses range from 240-500 MeV and are not much lighter than the physical kaon, which is 
integrated out in SU{2) xPT. Furthermore, the largest of the taste-splittings on the coarse 
lattices is not much smaller than the kaon mass [a^A/ ^ (460MeV)^], though on the fine 
lattices it is a factor of 2.7 times smaller [a^A/ ^ (280MeV)^]. 

There have been several studies showing that NLO SU{3) mixed-action and staggered 
chiral perturbation theory accurately describe the lattice discretization effects even at the 



171,1561,1571. For example, the 



rather large taste-splittings on the coarse MILC lattices 
agreement between the NLO mixed-action SU (3) xPT prediction of the scalar correlator at 
large times and the lattice data is excellent over the range of masses and lattice spacings 



used in the present calculation of Bk jl7|. The agreement is good to within the statistical 
accuracy of the scalar correlator data, which is about 10%. This agreement is a highly 
non-trivial test because the heavy taste-singlet r] plays the dominant role in modifying the 
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continuum xPT form of the correlator. The MILC collaboration also finds that the scalar 
correlator in the staggered theory is well-described by staggered chiral perturbation theory, 
and that when the low energy constants in the prediction for the scalar correlator are allowed 



16|. 



to vary, they agree with those determined in fits to the light pseudoscalar sector 

The statistical errors on are ~ 0.5% — 2% for most of our data points. It is now well- 
established that NLO xPT does not describe quantities such as pseudoscalar masses, decay 
constants, or to percent-level accuracy at the physical kaon mass, nor is it expected 
to based on power counting. Our data set confirms this picture for B^- In order to get 
good fits to our Bk data in the region of interest we must include next-to-next-to-leading 
order (NNLO) analytic terms. Fits without these terms give terrible correlated x^/d-o.f.'s 
and miniscule confidence levels. The two-loop NNLO logarithmic corrections, however, are 
not known for Bk. These expressions would also have to be modified to account for the 
staggered sea sector, though, given our experience with the one-loop modifications due to 
the mixed action, this is likely a small effect. In the region where the NNLO analytic terms 
that we have added are important, we expect the NNLO logarithms to vary slowly enough 
that their effect is well approximated by the analytic terms. We vary the assumptions we 
make for the (thus far unknown) NNLO behavior in order to estimate a systematic error for 
the chiral extrapolation. When we apply the same approach to our calculation of /tt, Jk, 



and their ratio 



231], we find systematic errors due to the chiral extrapolation of a similar size 



to other lattice determinations 26|] and in excellent agreement with phenomenology. This 
good agreement between our lattice calculations and known quantities lends confidence in 
our methods for quantities such as Bk that are not known from experiment. 

Our approach to chiral fits using SU (3) xPT plus higher-order analytic terms is related 
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in the 



to the approach of other groups using SU (2) xPT with the kaon integrated out 
following way. The chiral logarithms due to pion loops are common to SU{3) and SU(2) 
XPT, and are the dominant non-analytic contribution when the strange quark is much 
heavier than the two lightest quarks. In order to get acceptable fits, we need to introduce 
polynomial dependence in the valence quark mass at higher order than NLO in the chiral 
expansion, as discussed above. This is necessary to describe the data in the region where 
the SU{3) xPT is not especially convergent, and higher-order corrections are important. 
Nonetheless, our simulations interpolate about the physical strange quark mass, so the 
higher-order dependence of the SU (3) expansion, with terms involving kaon and eta masses. 
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is expected to be well described. Towards the physical value of Bk, where we extrapolate 
in the light quark masses, the heavier meson masses in the SU (3) expansion decouple and 
the SU{3) form becomes that of SU{2). In this decoupling region, the taste-breaking in the 
heavier mesons containing the strange quark reduces to terms analytic in and the light 
quark masses, which are included in the fits. Thus, we expect our fit function to describe 
our light-mass data well, and we expect the additional analytic terms to capture higher- 
order effects in the heavier-mass region where terms beyond NLO are necessary. Recent 
work by the MILC collaboration corroborates this approach. They find excellent agreement 
between SU (2) and SU{3) staggered chiral perturbation theory approaches for and light 



quark masses 59|, |60| . The vt-tt scattering results of NPLQCD also show good agreement for 
';he 1 = 2 scattering length a^^^ between SU{2) and SU{3) NLO MA^PT determinations 

There are six new continuum NNLO analytic terms for B^, as well as NNLO terms that 
modify the NLO constants C1-C4 by terms proportional to a^. We include only a subset of 
the NNLO terms that are needed to obtain good correlated values, and include the others 
in alternative fits for systematic error estimation. The number of new continuum low energy 



constants can be constrained using CPS symmetry [6l|, chiral symmetry, and the fact that 
there is only one mass scale in the tree- level diagrams with the external kaons at rest. The 
new continuum analytic NNLO contributions to B^ are 

4«)(2mi^ + ml^\ d y^-f^^\ 2ml^ + m^J. (62) 

We also test for higher order analytic terms proportional to a^. We find an improvement 
to the fit when including a term of the form a^m^^. Our systematic error estimate includes 
the effects of generic NNLO non-analytic terms on the extrapolation and an estimate of the 
size of NNLO taste- violations not accounted for in the fitting procedure. 

Figure [11] shows our preferred fit to the data using NLO partially quenched MA^PT 
supplemented by some of the above NNLO analytic terms. In order to obtain a fit with a good 
correlated confidence level, we include the NNLO continuum analytic terms proportional to 
di, d2, d^ and d^ in Eq. ( l62l) plus an NNLO analytic term containing discretization effects, 
da2m2y{8/ f^y)a^mly. We fix the following parameters in the fit: the tree-level (continuum) 
coefficients fi^w and /^stag, the decay constant f^y, and the taste-splitting a^Aj. We take for 
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FIG. 11: Bk versus degenerate light valence quark mass ri{mx+mj-es) on different ensembles. The 
fit lines are the same partially quenched contours as the data. Note that all of the coarse fit lines 
lie on top of each other, and likewise for the fine, indicating very little sea quark dependence. The 
band is the degenerate quark mass full QCD curve {rrix = niy = mi = rrih) in the continuum limit. 
The y-intercept of the full QCD curve gives the low-energy constant Bq, which is the value of Bk 
in the SU{2>) chiral limit. 




TABLE V: Tree-level //-parameters, defined in Eqs. (j46|) -(|47 p . on the coarse and fine MILC lattices. 
The staggered coefficient, ri/Xstagi is independent of the valence sector Q], while the domain- wall 
coefficient, ri/^dw) depends upon the choice of HYP-smearing [l^]. 



a(fm) 


ri/Xstag 




0.12 


6.234 


4.13 


0.09 


6.382 


3.83 
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FIG. 12: Comparison of higher order corrections for Bk- The rightmost point on the graph 

corresponds to ~ ms/2. 



the parameters /idw and /istag the values obtained from fits to the hght pseudo-scalar masses 
squared to the tree- level forms given in Eqs. fllBl) and (jlT]). This accounts for higher-order 
chiral corrections and is more accurate than using fi obtained in the chiral limit (which is 
found by fitting to the one- loop pseudoscalar mass and decay constant expressions), giving 
a better approximation to the pion mass squared at a given light quark mass. The values 
for /idw and //stag are given in Table |Vl We take the decay constant f^y, which appears as 
the inverse square in the coefficient of the chiral logarithms, to be the physical Jk = 156.5 
MeV 2J] for our preferred fit, though we vary fxy in order to estimate the systematic error. 
We use the value for the taste-singlet splitting a^Aj obtained by the MILC Collaboration in 
Ref. [3] and given in Table HVl Given these choices, our preferred combined chiral-continuum 
extrapolation fit function contains only ten free parameters. 

Figure [TT] shows only the degenerate valence mass points and the corresponding partially 
quenched fit lines, although the fit includes nondegenerate masses as well. The heaviest 
valence kaon masses included in this fit are slightly larger than the physical kaon mass. 
We restrict the degenerate valence "kaon" masses to below 500 MeV, but we allow slightly 
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heavier non-degenerate valence kaons up to masses of around 600 MeV in order to interpolate 
about the physical strange quark mass. In the sea sector, we restrict the taste-pseudoscalar 
pions to be less than 550 MeV on the coarse ensembles and less than 500 MeV on the fine 
ensembles. Our lightest degenerate valence "kaon" is ~ 230 MeV, while our lightest taste- 
pseudoscalar sea pion is ~ 240 MeV. Given these mass restrictions, the number of data 
points in our preferred fit is 69, which is more than sufficient to constrain ten parameters. 

Although most of our degenerate-mass data points are far from the physical kaon mass, 
including this data in the fit allows us to constrain the parameters of the SU{3) chiral La- 
grangian and to study the convergence of SU{3) xPT. The (cyan) band in Fig. [TTl shows the 
full QCD curve with statistical errors in the SU{3) {m^ = my = rrii = rrih) and continuum 
limits that is obtained from our preferred fit. In order to examine the convergence of SU{3) 
xPT, we plot the separate contributions to the degenerate SU{3) curve through LO, NLO, 
and "all orders" in Fig. [121 The right-most part of the x-axis corresponds to ~ 772^/2 (i.e., 
a 500 MeV pion), where we do not expect xPT to be especially convergent. Because we 
are interpolating in the quark mass in this region, we expect the "all orders" curve to be 
accurate, with the NNLO terms approximating the correct higher order behavior. Closer to 
the physical pion mass, however, the NLO contributions are the dominant corrections, and 
the particular choice of NNLO analytic terms has little impact on the fit result in this light 
quark-mass region as long as the fit has a good correlated confidence level. 

It should be noted that the relative contributions of NLO versus NNLO terms can change 
significantly, depending on whether one uses bare expansion parameters or the physical 



masses and decay constants in the one-loop expressions. Following Bijnens 62|, we prefer to 
use the physical masses and decay constants because the chiral logs are created by particles 
propagating with their physical momentum, and (though not relevant for Bk) thresholds 
appear in the right places at each order in perturbation theory. Using the bare parameters 
in our fits yields poor confidence levels, and we do not consider them further. Using the 
physical parameters in the expansion, we find the results of our fits to be consistent with the 
expectations from chiral power-counting. This is consistent with the findings of the JLQCD 
collaboration that using a resummed physical expansion parameter significantly improves 
the ability of NLO (or NNLO) xPT to describe data at heavier masses near the kaon by 



accelerating the convergence of the chiral expansion 



63|. The NPLQCD collaboration also 



uses the physical parameters in the SU (3) MA^PT expansion for their determination of the 
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vr-vr scattering length, and they find good a gree ment between their data and the leading 



order form, which is completely predicted 
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57|. When NPLQCD include the NLO SU{3) 



corrections, they also find good fits to their lattice data. Despite these successes of the 
SU (3) formalism, it would be valuable to continue our study with the complete NNLO 
formula once it is available. It should also be noted that the systematic error in the chiral 
extrapolation to the SU{3) chiral limit is fairly large, since the simulated strange quark 
masses in the sea are close to the physical value. Thus, there is a large systematic error in 
the value of the leading order term Bq, and the picture in Fig. [T2] may change appreciably 
once we can better constrain the approach to the SU{3) limit. Given the central value of 
our current best fit, however, the convergence of the chiral expansion appears reasonable. 

We obtain a value of i?^'^(2GeV) = 0.5273(64) from our preferred fit when the matching 
factor is calculated using NPR, where the error is statistical only. We take this result as our 
central value. For comparison, the value for B^ obtained using lattice perturbation theory 
to compute the matching factor is Bf^'^{2GeY) = 0.541(6). The result of our preferred fit 
is shown in Fig. [131 All data points used in the fit are shown. The upper band is the full 
QCD continuum extrapolated curve with the strange quark fixed to its tuned value. The 
lower band is the degenerate quark mass, fuU-QCD band, as in Fig. [TTl The extrapolated 
value of Bk at the physical quark masses with statistical errors is shown as an "X" with 
solid black error bars. 

Table IVTl shows the low-energy constants determined in our preferred fit, with statistical 
errors only. We do not attempt to estimate a systematic uncertainty in these parameters 
because this is not necessary for determining Bk at the physical quark masses. We expect 
that the systematic uncertainties, however, will be large given the size of the extrapolation 
to the SU (3) chiral limit. We present the values only to illustrate a few important points. 
Discretization errors in our data are small; this can be seen from the size of the parameters 
Ci and daimly- Further, we do not observe any clear sea-quark mass dependence in our data, 
as shown by the fact that C4, d^, and dg are zero within errors. We estimate the systematic 
uncertainty in Bk due to the choice of chiral and continuum extrapolation fit function in 
the following section. 
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FIG. 13: Bk versus light valence quark mass ri{mx + Tn^cs)- All data points used in the fit (same 
fit as in Fig. [TT]) are shown. The upper band is the full QCD continuum extrapolated curve with 
the strange quark fixed to its tuned value. The lower band is the degenerate quark mass, full-QCD 
band, as in Fig. [TTJ The extrapolated value of Bk at the physical quark masses with statistical 
errors is given by the "X" with solid black error bars. 

V. SYSTEMATIC UNCERTAINTIES 

In the following subsections, we examine the uncertainties in our calculation due to 
the chiral/continuum extrapolation, scale and light quark mass uncertainties, finite volume 
effects, and uncertainties in the matching factor Zbi^- 

A. Chiral and continuum extrapolation errors 

We estimate the systematic error in the chiral extrapolation by varying the fit function 
used to extrapolate the data over a variety of different reasonable choices and taking the 
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TABLE VI: Fit parameters obtained in our preferred fit. Errors shown are statistical only, and do 
not include the extrapolation uncertainty. The coefficient Bq in the top panel is the only leading- 
order low-energy constant, the coefficients in the middle panel are the NLO low-energy constants, 
and the coefficients in the bottom panel are the NNLO parameters included in the fit. Definitions 
of the parameters are given in Eqs. (I49p and ()62p . 



Bo 


0.338(12) 


Cl 


0.057(25) 


C2 


0.00531(60) 


C3 


0.00033(15) 


C4 


0.00003(15) 


di 


0.351(23) 


d2 


-0.004(38) 




-0.006(27) 


de 


0.006(35) 


da^mly 


-0.00049(29) 



spread between them. By reasonable, we mean theoretically motivated fits that also describe 
the data with good confidence levels, determined by the correlated per degree of freedom. 
These fits always involve the known one-loop mixed action chiral logarithms, since including 
them incorporates the leading non-analytic dependence on the light quark masses. We 
consider variations of the fit function by including different types of terms beyond NLO 
that still give acceptable confidence levels and in order to determine the systematic error in 
the chiral extrapolation. We also vary other assumptions, such as the values of parameters 
used in the prediction for the chiral logarithms, and the lattice spacing dependence of the 
continuum extrapolation. Each of these variations is addressed in turn. 

We combine the chiral and continuum extrapolations using MA^PT to control the ap- 
proach to physical light quark masses and to the continuum. Combining the data sets on 
coarse and fine lattices, we have seven different valence quark masses and nine different sea 
quark mass combinations. Our valence kaons range from around 600 MeV down to as light 
as 230 MeV. These lighter kaons are useful for constraining the low energy constants of 
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the chiral Lagrangian. As can be seen in Fig. [121 even at the physical kaon mass near the 
rightmost part of the plot, the chiral behavior is mostly accounted for by a combination of 
leading order and NLO terms. The combinations of sea quark masses include values of the 
simulated strange sea quark above and below the physical strange quark mass on both the 
coarse and fine lattices, allowing us to interpolate in the strange sea quark mass. 

The light sea quark masses used in our simulation are as low as m^/lO, and this translates 
into a taste-Goldstone pion (the lightest of the staggered pions) of around 240 MeV. In 
the chiral extrapolation of B^, however, the only sea pion mass that appears in the NLO 
MA^PT expression is the taste-singlet pion mass (the heaviest of the staggered pions), the 
lightest of which in our simulations is still a rather heavy 370 MeV. Fortunately, the sea 
quark contribution to the NLO chiral logarithms vanishes for degenerate valence quarks, 
and gives only a small contribution for nondegenerate valence quarks (the region of interest 
for the physical kaon). This is in part because the terms that contain the taste-singlet pion 
are suppressed by a factor of l/Ngca- The taste-breaking is corrected for explicitly through 
NLO, and partially at NNLO by including analytic terms proportional to and a^rn^y. 

As described in the previous subsection, we take as our central value the result of a fit that 
includes all of the terms through NLO, the NNLO continuum analytic terms proportional 
to c/i, c?5 and in Eq. (l62l) . and an NNLO analytic term containing discretization effects 
proportional to a'^m'^y. We estimate the systematic uncertainty due to the chiral extrapo- 
lation by including additional NNLO analytic terms and taking the difference between the 
new value of Bk and that obtained with the preferred fit. The largest contribution to the 
systematic error comes from adding terms quadratic in the sea quark masses to the above 
fit and taking the spread between the two. Although we analyze ensembles with different 
light and strange sea quark masses, we do not observe any clear sea-quark mass dependence 
in our data. When we include the terms proportional to ^3 and in Eq. (!62|) to our fit the 
result shifts to 0.5175(72), yielding a difference of 1.9%. We take this as the error due to 
approximating higher order terms in the chiral expansion. 

We also consider NNLO non-analytic terms of the generic form log(m^) to estimate 
higher order effects. A term of the form m^[log(m^/A^)] appears at NNLO [in both SU{2) 
and SU{3) xPT], but it is subleading in the chiral expansion and should therefore have a 
smaller impact than the terms that we are already including as rrix approaches the physical 
m^r. In order to test that the effects of such a term are indeed small, we added this term to 
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our preferred fit leaving its coefficient as a free parameter. We obtain a small coefficient for 
such a term, which leads to a slight 0.9% shift upwards in our central value. This is within 
our estimate of the error due to approximating higher order terms in the chiral expansion. 

Higher-order taste-breaking is considered as well. If we set the staggered singlet taste- 
splitting to zero in the NLO logarithms this amounts to using the continuum-like expression 
appropriate to a purely domain- wall simulation; if we do this it shifts our chiral extrapolation 
to Bk at the physical quark masses by only 0.7%. Thus, the discretization effects particular 
to the use of staggered quarks in the sea sector are small. The formulas we use do not 
explicitly remove taste-breaking at NNLO. Note, however, that the higher order analytic 
term a?m?^y is needed to get good confidence levels in our fits, and is expected to absorb 
some of the higher-order taste- violating effects, especially those in the kaon sector of SU (3) 
xPT in the limit that the kaon can be integrated out. We estimate the residual effect of the 
higher-order taste-breaking by varying the splittings used for the sea-mesons in the analytic 
terms [Eq. ( l62l) ] over the full range of staggered meson masses. When we do this, we find that 
the central value decreases by 1.8%, which is within our systematic error due to neglecting 
higher order terms in the chiral expansion. 

We consider other variations to the fit, but they lead to a much smaller shift in the 
central value. Although we fix the tree- level (continuum) coefficients /idw and /istag? the 
decay constant f^y, and the taste-splitting a^A/ in the NLO chiral extrapolation formula, 
we vary them within their statistical uncertainties in order to estimate their contribution to 
the error in Bk- The staggered and domain- wall //tree parameters are well-determined, and 
their error is negligible in Bk- In the chiral fit used for our central value we take the decay 
constant f^y-, which appears as the inverse square in the coefficient of the chiral logarithms, 
to be the physical /x- We vary this coefficient between and /x as an additional way 
of estimating higher order corrections. Note that a change in the coefficient of the chiral 
logarithms will change the other fit parameters, so that this produces a much smaller effect 
than simply changing the overall coefficient of the chiral logarithms in the final result by a 
factor of (/a-Z/tt)^ ~ 1.4. The change of fxy from fK to f-,, leads to a 0.1% shift in the value 
for Bk- The approximately 10% statistical error in a^A/ leads to a similarly negligible error 
in Bk- 

Although we include terms proportional to in the preferred fit, there is some ambiguity 
(with only two lattice spacings) in the dominant source of discretization errors, which may 
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FIG. 14: Same fit as Fig. [T3l but witli a subset of tfie data points for illustration. In this case, the 
valence masses are nondegenerate, with the heavier mass fixed close to the strange quark mass. 
The light-quark mass is the lightest simulated on each ensemble. The fit curve is the full QCD 
continuum extrapolated curve with the strange quark fixed to its physical value. The extrapolated 
value of Bk is shown, including the statistical error (solid error bar with X) and the systematic 
error due to the chiral extrapolation, combined with the statistical error in quadrature (dashed 
error bar). The dotted error bar (star, slightly offset) shows the total error for Bk- 

be purely corrections, taste-breaking terms proportional to a^a^, or chiral symmetry 
breaking terms proportional to mrcsO^- It is also possible that the different sources lead 
to discretization effects of the same size. We therefore vary the change in the effective 
between coarse and fine lattice spacings, taking the resulting spread in the value of as 
part of the systematic error. If discretization effects decrease as mresO^ or as a^a^ then they 
should go down by about a factor of three from our coarse to fine lattices. If they decrease 
as they should decrease by about a factor of two. This difference leads to a 0.3% change 
in the continuum extrapolated central value. 

In summary, the largest source of chiral-continuum extrapolation error comes from un- 
certainty in the sea-quark mass dependence. The parametric uncertainty on quantities used 
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as inputs in the chiral logarithms is neghgible, and the residual errors due to the lattice 
spacing dependence in the continuum extrapolation are small. Adding the uncertainty due 
to approximating higher order terms in the chiral expansion and the residual continuum 
extrapolation error in quadrature, we quote a total systematic error due to the chiral and 
continuum extrapolation of 1.9%. 

Figure [14] illustrates the chiral extrapolation error. This is the same fit as that shown 
in Fig. [131 but with a subset of the data points. In this case the valence masses are non- 
degenerate, with the heavier mass fixed close to the strange quark mass. The light quark 
mass is the lightest simulated on each ensemble. The fit curve is the full QCD continuum 
extrapolated curve with the strange quark fixed to its tuned value. We extrapolate the light 
valence quark mass to the physical d quark mass, while we extrapolate the light sea quark 
mass to the average of the u and d quark masses. The band shows the full QCD curve ending 
at the full QCD d quark mass. The error bar is centered on the final result, which has a 
small (not visible) shift due to setting the light sea quark mass equal to the isospin-averaged 
quark mass. The extrapolated value of Bk is shown, including the statistical error (solid 
error bar with X) and the systematic error due to the chiral extrapolation, combined with 
the statistical error in quadrature (dashed error bar). The dotted error bar (star, slightly 
offset) shows the total error for Bk including the matching error. 



B. Scale and quark mass uncertainties 



In order to convert lattice quantities into physical units we use the MILC Collabora- 
tion's determination of the scale, ri, where ri is related to the force between static quarks. 



r\F{ri) = 1.0 33|, |3J]. The ratio ri/a can be calculated precisely on each ensemble from 



the static quark potential. We use the mass-independent prescription for ri described in 
Ref. 2J] so that all of the mass dependence is explicit in MAxPT and none is hidden in the 
scale-fixing scheme. In order to fix the absolute lattice scale, one must compute a phy sical 
quantity that can be compared directly to ex per iment; we use the T 2S-1S splitting 35j and 



the most recent MILC determination of /^r 



24| . The combination of the T mass-splitting 



and the continuum-extrapolated ri value at physical quark masses leads to the determina- 
tion r?''^' = 0.318(7) fm [Sel. The use of ^ to set the scale yields rf^' = 0.3108(15)(^?6) 
fm 2^. This difference between the two scale determinations leads to a systematic error in 
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our result for Bk- Since Bk is a dimensionless quantity, the scale enters only through the 
quark mass determinations. We determine the light valence quark masses using MA^PT fits 



to the pseudoscalar masses and decay constants, as described in Ref. [23|. The sea-quark 



masses are taken from the most recent update of the MILC pseudoscalar analysis 
MILC finds for the bare staggered quark masses in ri units 



24 



64| 



24|. 



^^^stag ^ ^q3 ^ 3.78(16), (63) 

nmf^s X 10^ = 102(4) , (64) 

where rh = {rriu + mrf)/2. Although the masses are in dimensionless ri units, they are scale 
and scheme dependent quantities. The scheme, of course, is the improved staggered lattice 
action used in the MILC simulations. The scale is the fine lattice scale of ~ 2.3 GeV, but 
with discretization effects removed by fits to multiple lattice spacings using rooted staggered 
xPT. Both our fits and the MILC fits use as inputs from experiment the averaged meson 
masses with electromagnetic effects removed as well as possible. We (and MILC) take for 
the squared meson masses m| and m|>. 



ml = mlo, (65) 
= \i^KO + - (1 + AE){ml+ - m^o)), 

where 1 parametrizes corrections to Dashen's theorem. 

We find from our mixed action chiral fits to the pseudoscalar sector the values for the 
bare domain-wall quark masses (also evaluated at the fine lattice scale) 

nm X 10^ = 5.87(8)(41), (66) 
nm, X 10^ = 168(2)(8), (67) 

where the first error is statistical and the second is systematic. Following MILC, we also 
obtain the masses of the two lightest quarks. Given m^, we can obtain m„ by extrapolating 
not to the mass of the K but to the mass of the (with EM effects removed). We take 

(m^+)QCD = m'^K+ - (1 + SE){ml+ - mlo), (68) 

where 5^; = 1, which corresponds to vanishing EM corrections to the mass. We then 
obtain 

rim„ X 10^ = 3.7(22)(7), (69) 

nnid X 10^ = 8.0(27)(7), (70) 
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where again the first error is statistical and the second is systematic. 

It is useful to observe for our Bk error analysis that the systematic error for the domain- 
wall rris is dominated by the scale error. However, the error in (needed since a kaon 
is an sd state) is dominated by statistical uncertainty. Thus we can treat the errors from 
the s and d quark masses as uncorrelated. Note that all of the above masses are the bare 
lattice masses, so no error has been included for the renormalization needed to match to 
a continuum scheme like MS. The bare quark masses are sufficient for the purpose of 
calculating B^- The error in rim^ leads to an 0.8% uncertainty in B^, while the error in 
rirud leads to an 0.2% error in Bk- The errors in the sea quark masses produce a negligible 
uncertainty in B^- Combining these errors in quadrature we obtain an error due to scale 
and quark mass uncertainties for B^ of 0.8%. 



C. Finite volume error 

The finite volume error is estimated using one-loop finite volume MA^PT joj. We have 
simulated at fairly large volumes, such that rriT^L > 3.5, and we have corrected our data 
using the appropriate one-loop MA^PT expressions, which are never larger than 0.6%. 
There could still be non-negligible residual finite volume corrections, however, as numerical 
studies by the MILC collaboration of and show that the one-loop xPT corrections 
can be off by as much as 50% for similar simulation parameters using staggered quarks 24 1. 
Even so, given that the largest finite-volume correction to any individual data point in our 
analysis is 0.6%, we expect the residual corrections to be only as much as 0.3%. However, 
in order to be conservative, we take the entire 0.6% as our total finite volume error. 



D. Renormalization factor uncertainty 

In this section we estimate the systematic uncertainty in B^ due to the nonperturbative 
determination of the renormalization factor Zb^- We consider several sources, discussing 
each in turn. 
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1. Chiral extrapolation fit ansatz 



In order to remove explicit chiral symmetry breaking contributions to Zbj^ from operators 
such as those in Eq. (EZ]), we first extrapolate to the chiral limit at fixed values 

of (ap)"^. Although we choose to use a fit function that is linear in the light sea quark 
mass and quadratic in the valence quark mass, we can obtain an equally good correlated 
confidence level using a fit function with even more terms. We must therefore consider the 
systematic uncertainty introduced by the choice of chiral extrapolation fit ansatz. We do so 
by adding a quadratic term in the light sea quark mass to the fit function and re-doing the 
chiral extrapolation at each value of (ap)^. We then re-compute Z^^(2 GeV) and take the 
difference between this result and the central value to be the systematic error. This leads 
to an uncertainty in 2'^^'™'^'''^'^ (2 GeV) of 0.0063, or ~ 0.7%, on the coarse lattice and an 
uncertainty of 0.0111, or ~ 1.2%, on the fine lattice. The addition of yet another term cubic 
in the valence quark mass produces a negligible difference in Zbj^- We take the larger, 1.2% 
difference, to be the uncertainty in Zbj^ from the choice of chiral extrapolation fit function. 



2. Strange sea-quark mass dependence 

When we extrapolate 2'^^^°^ to the chiral limit at fixed values of {apY, we do not, in 
fact, take the value of the strange sea quark mass to zero. This is because, within statistical 
errors, is independent of the strange sea quark mass. We can explicitly calculate the 
strange sea-quark mass dependence, however, by taking the chiral limit of z^j^^'^^ at fixed 
{apY using a function that is linear in the sum of the sea quark masses and quadratic in the 
valence quark mass on the coarse lattices, so we make the replacement 2mi — )■ (2m; + mh) in 
Eq. (gS]). We find that this leads to a difference in z^s,coarscj.2 q^^) of 0.0027, or ~ 0.3%, 
and take this to be the uncertainty in due to the nonzero strange sea quark mass. 



3. Aa - Ay 7^ 

The use of exceptional kinematics in our nonperturbative renormalization factor calcu- 
lation leads to a difference between Ka and Ay of ~ 1% at nonzero quark masses and 
p ^ 2 GeV. This is shown in Figs. H] and [51 Because we do not know a priori which of the 
two quantities has less contamination from chiral symmetry breaking, we use the average 
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{Aa + A\/)/2 to determine the central value for Zbjc In order to estimate the systematic 
uncertainty that is introduced by this choice, we also calculate using for the normal- 
ization. This leads to a difference in ^^^^^"'^'^^^ ^2 GeV) of 0.0084, or ~ 0.9%, on the coarse 
lattice and a difference of 0.0112, or ~ 1.2%, on the fine lattice. We take the larger, 1.2% 
difference, to be the uncertainty in Z^^. due to chiral symmetry breaking between and 
Av. 



4- Mixing with wrong- chirality operators 

The use of exceptional kinematics also leads to mixing between the standard model 
operator Ok, which has a VV + AA chiral structure, and other operators of different chi- 
ralities that do not contribute to — mixing in the standard model. Although the 
size of the mixing coefficents, shown in Figs. [TSHTSl are small, the matrix elements for 
the wrong-chirality operators diverge in the chiral limit and are much larger than the de- 



sired matrix element 
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Thus a small mixing coefficient can still potentially lead to a 
non- negligible error in Zbi^- Fortunately, we can bypass this concern by computing the 
mixing coefficients at non-exceptional kinematics. Theoretically, we expect their size to be 



of 0{{amresY) ~ 10"^ |49|, l52|. Numerically, we find that all of the mixing coefficients are 
consistent with zero on both the coarse and fine lattices, as shown in Figs. [T91 - I261 Because 
the contribution to Bk in the MS scheme from each wrong-chirality lattice operator is in- 
dependent of the lattice scheme initially used to obtain the mixing coefficients, we conclude 
it is safe to neglect them in our calculation of ^b^, despite the fact that we are using ex- 
ceptional kinematics. We therefore do not add any systematic uncertainty to Zbj^ due to 
operator-mixing. 



5. Perturbative matching and running 

Although we compute Zb^^ in the RI/MOM scheme nonperturbatively, we must still 
convert its value to the SI scheme to remove lattice discretization effects and ultimately to 
the MS scheme using 1-loop continuum perturbation theory. This introduces uncertainty 
into Bk due to the omission of higher-order terms. Because the true truncation error cannot 
be known without the computation of the next term in the perturbative series, we consider 
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several ways to estimate the uncertainty here. 

The first is to muhiply the largest individual 1-loop conversion that is used by an ad- 
ditional factor of a^^. We determine from data in the momentum window 2 GeV < 
p < 2.5 GeV; thus the largest value of a^^(p) used is that at 2 GeV. The largest correction 
comes from the conversion between the RI/MOM scheme and the SI scheme, and leads to 
the following estimate of the truncation error: 

^MS(2 GeV) X ^^^^^j^JS/mom = 0-0188, (71) 

or ~ 2%. 

The second is to take the size of the entire 1-loop correction from the RI/MOM scheme 
to the SI scheme to the MS scheme: 



afS(2 GeV) 

4^ RI/MOM ^MS 



USlMOM - JS.) = 0.0204, (72) 



which also leads to an estimate of ~ 2%. Because the conversion factors, however, are only 
known to 1-loop, they are particularly sensitive to the scale at which they are evaluated. 
We did not attempt to determine an optimal scale for the process, for example using the 
BLM prescription 4^, and must therefore estimate the error due to scale ambiguity. The 
standard, although somewhat arbitrary, prescription used in the continuum literature is to 
take the variation in the quantity when the scale fi is varied between 2/i and /i/2. For our 
case, this leads to the estimates 

^5r^'^-/«OM-^S)= 0.0152. (73) 

^^^(jShMOM - ^S) = 0.0327. (74) 

Thus, the 1-loop correction can be as large as ~ 3%, if we use a scale of 1 GeV. 

The third is to take the difference between Zbj^ determined using the nonperturbative 
Rome-Southampton approach and using lattice perturbation theory. Each method for com- 
puting Zbj^ relies on 1-loop perturbation theory, but involves a different series expansion, so 
one does not know a priori which leads to a faster converging series and smaller truncation 
error. Thus having two independent calculations of Zbj^- provides a valuable independent 
crosscheck. To estimate the error, we replace the values of Zbi^ determined using nonpertur- 
bative renormalization with those from lattice perturbation theory and repeat the extrapo- 
lation to the physical quark masses and the continuum. We obtain 5^^^(2GeV) = 0.541(6), 
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where the error is statistical only. We then take the difference between the resulting Bk 
and our preferred central value: 

\B^.^^ - Bjf^\/B^/^ = 0.027, (75) 

which is ~ 3%. This is comparable to, but slightly smaller than, the estimate from the scale 
ambiguity. Nevertheless, we think that it is a more reasonable estimate of the uncertainty, 
given that it comes from two independent perturbative computations. 

Finally, we consider the possibility that the scale dependence observed in (see Figs. M 
and [To]) is not due solely to discretization errors. Although we made this assumption when 
obtaining the central value for B^, some of the slope in versus {apY may in fact be 
due to the lack of higher-order terms in the matching factor, which is currently known to 
only 1-loop in perturbation theory. Our fourth method for estimating the truncation error 
is therefore to take the difference between Z^^ at 2 GeV and that obtained in the limit 
p — )■ 0. To estimate the resulting error in B^, we replace the values of used in our 
preferred analysis with those obtained without extrapolating Z^^^ to zero momentum, and 
repeat the combined chiral and continuum extrapolation. We obtain i?^^'^(2GeV) = 0.542(7) 
(statistical error only), which is remarkably close to the lattice perturbation theory value. 
We then take the difference between the resulting B^ and our central value: 

|5^^o - S^=2G-^|/5^^-^o = 0.028, (76) 

which is ~ 3%, and similar to our previous estimate. Although taking Z^^ directly at 2 
GeV provides a sensible alternate method for obtaining the renormalization factor, it still 
does not eliminate truncation errors from the lack of higher-order perturbative matching. 
Nevertheless, because the values of B^ obtained using this alternate nonperturbative de- 
termination and from lattice perturbation theory are so similar, we estimate the residual 
truncation error to be small: 

|5^=2GeV _ 5iPT|/^^=2GeV ^ q_qq2, (77) 

which leads to a negligible increase in the systematic error when added in quadrature. 
Since this method for estimating the perturbative truncation error leads to a slightly more 
conservative estimate than the third approach, we take 2.8% to be the uncertainty in 
due to the use of 1-loop perturbation theory. 
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TABLE VII: Error contributions to Bx from the nonperturbative renormalization procedure. Each 
source of uncertainty is discussed in Sec. IV Dl and is given as a percentage of Bk- 



uncertainty 



statistics 


0.7% 


chiral extrapolation fit function 


1.2 


% 


strange quark mass dependence 


0.3 


% 


chiral symmetry breaking 


1.2 


% 


perturbation theory 


2.8 


% 



total 3.4% 



6. Total uncertainty in Zbj^ 

We summarize the contributions to the "renormalization factor" uncertainty in Bk in 
Table IVIII and add them in quadrature. The ~ 2.8% error due to the use of perturbation 
theory is the largest single contribution to the total error in B^, and can only be reduced 
by a calculation of the necessary matching factors at 2-loops. 

VI. RESULT AND CONCLUSIONS 

We obtain the following result for Bk in the MS scheme at 2 GeV: 

B^{2 GeV) = 0.527(6)(10)(4)(3)(18), (78) 

where the errors are from statistics, the chiral-continuum extrapolation, scale and quark mass 
uncertainties, finite volume errors, and the renormalization factor uncertainty, respectively. 
The total error is ~ 4%, and the error budget is presented in Table IVIIIi It is often more 
convenient to use the scale-invariant parameter B^ in new physics analyses, for which we 
find the value 

= 0.724(8)(29). (79) 

Our 2+1 flavor lattice QCD calculation of Bk is the first to have all lattice sources of 
systematic uncertainty under control. The largest errors in our result for Bk come from the 
chiral-continuum extrapolation (1.9%) and from the determination of the renormalization 
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TABLE VIII: Total error budget for B^- Each source of uncertainty is discussed in Sec. |Vl and is 
given as a percentage of Bk- 



uncertainty 


Bk 


statistics 


1.2% 


chiral & continuum extrapolation 


1.9% 


scale and quark mass uncertainties 


0.8% 


finite volume errors 


0.6% 


renormalization factor 


3.4% 


total 


4.2% 



factor (3.4%). The former uncertainty can be improved by the addition of statistics and the 
use of more lattice spacings. The MILC collaboration has generated ensembles with a lattice 
spacing of a 0.06 fm which we plan to analyze in the near future. The latter uncertainty 
can be improved in several ways. The use of Landau gauge- fixed momentum-source prop- 



agators will reduce the size of the statistical errors in Zb^ |65|, |66|, and may consequently 
better constrain the extrapolation to the chiral limit. The use of non-exceptional kinematics 
will reduce the contamination from chiral-symmetry breaking 53| and also provide an alter- 
native nonperturbative renormalization scheme with independent truncation errors from the 
standard RI/MOM scheme 67|]. A calculation of the 2- loop continuum perturbation theory 
formulae needed to match Zbi^ in the RI/MOM scheme to in the MS scheme would 
allow for a better estimate of the perturbative truncation error of Zbj^ in the RI/MOM 
scheme. Nevertheless, our calculation of the matching factor in mean-field improved 
lattice perturbation theory provides a robust alternative to our nonperturbative determina- 
tion in the RI/MOM scheme since the systematic uncertainties are uncorrelated between the 
two methods. In particular, the difference between the two results allows for a more reliable 
estimate of the matching error than from the RI/MOM scheme alone. This is important 
because some errors, such as perturbative truncation errors, are difficult to estimate within 
a single scheme. The error in Bk from all sources except the renormalization error is only 
2.5% because the use of domain-wall valence quarks and staggered sea quarks allows us to 
control the remaining sources of uncertainty quite well. Thus, if the use of non-exceptional 



53 



kinematics or 2-loop continuum perturbation theory does reduce the matching error, we can 
obtain an even more precise determination of Bk without making any further improvements 
to the lattice calculation. 

Our result is consistent with the determination by the RBC and UKQCD Collaborations 
using 2+1 flavors of domain-wall fermions, = 0.720(13) (37) |3i], but our result has a 
smaller total error. The largest error in the RBC/UKQCD calculation is the 4% scaling 
uncertainty due to the use of only a single lattice spacing, which we reduce by using two 
lattice spacings. Our result also has smaller statistical errors because of the large number 
of available staggered gauge configurations. Our result has comparable matching errors to 
RBC/UKQCD because the dominant error in both calculations of Zbj^^ is from the use of 
precisely the same one-loop continuum perturbation theory results when converting from 
the RI-MOM scheme to the MS scheme. Our error estimate is slightly more conservative, 
however, because we take the difference between the renormalization factors determined 
using NPR and using lattice perturbation theory to be the error due to the omission of 
higher-order terms. 

Our result is 1.9-cr lower than the value currently preferred by the global unitarity triangle 



analysis, Bk^ 0.92 ± 0.10 
Soni in Ref. 
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which comes from an update of the work of Lunghi and 



69| using the latest determinations of all of the input parameters. The tension 



with the standard model is enhanced by the inclusion of the correction factor derived by 
Buras and Guadagnoli 70|, iTlJ, which raises the location of the ex band. The uncertainty 
in the standard model constraint on Bk is ~11% . This is largely due to the error in the 
CKM matrix element |Vcfe|, which is known to ~ 2% accuracy, but enters the constraint from 
Bk on the unitarity triangle as the fourth power. Thus the error in \Vcb\ must be reduced 
in order to maximize the constraint on new physics from neutral kaon mixing. Fortunately 
work on improving the exclusive determination of \Vcb\ is ongoing by the Fermilab Lattice 
and MILC collaborations 72[, and work on improving the inclusive determination of \Vcb\ is 
in progress by Becher and Lunghi 73 1. 

Lattice QCD calculations of the hadronic weak matrix element Bk that incorporate the 
effects of the dynamical up, down, and strange quarks can now reliably control all sources of 
uncertainty. Because our result for Bk is consistent with the determination of the RBC and 
UKQCD Collaborations, one can safely average the two values (taking correlations between 
systematic errors into account) for use in future unitarity triangle analyses. There is already 
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a hint of the presence of new physics in the quark flavor sector as indicated by the tension 
between the unitarity triangle constraints from ex and sin(2/3) 69j. We expect the errors 
in both lattice QCD calculations of to be reduced in the future, such that indirect 
CP- violation in the kaon system will play a valuable role in the search for new physics. 
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Appendix A: Feynman rules for lattice perturbation theory 

In this appendix we present the Feynman rules and the integrals needed to calculate 
Zbj^ to one-loop in lattice perturbation theory with a Symanzik-improved gauge action and 
HYP-smeared domain wall quarks. 
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1. Gluon propagator 



The gluon propagator for the Symanzik-improved gauge action used by the MILC Col- 
laboration is 

kii.kij 



2\-2 



I -a) 



2\2 ■ 



where 



A,.{k) = A,^{k) = {l-S,,)A{kr' 



{ky - 



c^k'UY.ki + k-Y.'^ki\ 

\ p p¥=f^,'^ / 

+ci f(E^p')' +^'E^p' J2^r+ (kr n ^2) 



+^cil{Pr+2j2kt-k'Y.k'r 



P T + P 



with ci 



' \2ul ' 



(mo is the fourth root of the plaquette) and 
fc, = 2sin^,^2 = 5^^J. 



(Al) 



(A2) 



(A3) 



(A4) 



Without loss of generality, we adopt the Feynman gauge a = \. The above propagator is 
that of the tree-level (tadpole) improved gauge action 7J, |75| . The gluon propagator in the 
improved case is significantly more complicated then that from the Wilson plaquette gauge 
action, where ^pW™**" _ _ ^\yQ action used in the generation of the MILC ensembles 
is further improved through 1-loop, but this additional improvement introduces corrections 
of higher order than 1-loop in 2'^^, and is not needed here. 



2. Domain wall propagator 



For the domain- wall propagator, we make use of the results of Ref . 43|] . There are three 
types of domain-wall quark propagators. The first connects general flavor indices: 



N 



S{v) 



St 



Y^{-ilp sinp^ + W- + mM~)su Gr{u, t)PR 

u=l 

N 

+ sinp^ + W+ + mM+)su Gl{u, t)PL, (A5) 



u=l 
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where = (1 =t 75)/2 are projection matrices, s, t, and u are flavor indices, the mass 
matrices are 

f -W 1 \ ( -W \ 

-W ... 1 -W 



w- 



1 

-W J 



w- 



1 -w j 



(A6) 



M- 



(A7) 



and G_RL are 



A 



GR{s,t) = -[-(l-m2)(l-Tye-")e°(-2^+^+*)-(l-m2)(l-W^e")e- 



■a{s+t) 



-21^ sinh(a)(e"(-^+^"*) + e"(-^-^+*))] + Ae""!^"*! 



A 



(A8) 



6-^(5, t) = -[-(1 - m2)(l - yt/e°)e°(-2^+^+*-') - (1 - m')(l - iye-")e"(-''-*+2) 



-21^ sinh(a)(e"(-^+'-*) + e°(-^-^+*))] + Ae""!' 
cosh (a) = 



2W 



A 



2iysinh(a)' 
F = l- e'^W - m^{l - We-"), 

ly = I-M5 + 5^(1 -cos p^). 



(A9) 
(AlO) 

(All) 
(A12) 
(A13) 



In these formulas m is the domain- wall quark mass, and M5 is the do main- wall height. 
is the number of sites in the flfth dimension, i.e. the number of generalized flavors. 

The second propagator connects the physical quark fleld q with the fermion fleld of general 
flavor index, 

= 1(^7^ sin p^-m(l-iye-°))(e-"(^"'^)p« + e-(^-i)Pz,) 

+ ^[m(z7^ sinp^ - m(l - l^e"")) - F]e-'' {e-''^'-^^ Pr + c-^^^-'^Pl), 



(A14) 



(^(P, s)q{-p)) 



le -(^-^)p^ + e-"(^-i)Pfi)(z7Msinp^ - m(l - We-'')) 



^ /g-a(.-i)p^ + e-"(^-^)p^)e-"[m(z7^sinp^ - m(l - lye"")) - F]. 
F 

(A15) 
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The third propagator is that of the physical quark field 

S,{p) = {q{pH-p)) = (A16) 



which reduces in the continuum limit to 



ip+ [1 — WQ)m 
where Wq = 1 — M5. 

The following form of the propagators, where we perform the sum over generalized flavor 
indices, is useful for evaluating the vertex diagrams needed to renormalize 43| . 

SL,{p)^Y.L{s){^P{p,s)q{-p)) = ( ^(^ _ ^^^-a J (^^T^smp,. - (1 - We^) , (A18) 

00 

S.nip) ^ J2^q{pm-P,s))R{s) = Sl,{p), (A19) 

s=l 

SM = J2(^{pm-P^'))^^') = I _ ^^e-a sinp;. - m(l - VTe-")) , (A20) 

00 

Sr,{p) ^ = S,l{p) , (A21) 



N 



s=l 



with 

L{s) = W-'^Pn + ^/^^'^Pl), (A22) 
R{s) = {wf^Pn + wi'^-'^PL), (A23) 

where in the rightmost expressions we take the limit that the number of lattice sites in the 
fifth dimension is infinite. In principle, this limit should be taken after the momentum 
integral, but there is no difficulty with taking the limit first. We use the mean-field improved 
value = M5 — 4(1 — uq) throughout the perturbative calculation, as discussed in 

subsection IIII A[ 

3. Quark gluon vertices 



The quark gluon interaction vertices are 43 1 



V,%k,p),t = V,lik,p)Sst = -igT- (l,Vi,{k,p) + Vr,{k,p)) 5,,, (A24) 
V2tAk,v)st = <.(fc,p)5.* = ^/^{T",n(7A(A:,p)+Fi,(A:,p))5^,5,„ (A25) 
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where g is the couphng constant, T"" are the S'?7(3) generators, and 

Fi^(/c,p) = cos^(-A;^+p^), (A26) 
Vi^{k,p) = zsin^i-k^+p^). (A27) 

To account for the HYP-smearing of the valence quarks to the order we are working, 
the vertices must be modified by a form factor /i^^,. Since all gluons begin and end on 
fermion lines, the gluon propagator gets replaced by a more complicated propagator D^i^ — )• 



h^xDx^hua- The form factor is 



76| 



where 



h^x = S^xDx + (1 - S^x)G^x, (A28) 



Dx = l-d,J2sl + d2Yl - dstspl (A29) 

v,p^X 

G^x = s^,sxG^xik), (A30) 

G^xik) = d,- 4^4^ + ^3^, (A31) 



and = sin-^. In Eqs. flA29p - flA3ip . the indices /i. A, p and a are all different. The 
coefficients di are defined by 

2 4 
c?i = -ai(l + 02(1 + tts)), = -"102(1 + 203), ds = 8aia2a3, (A32) 

where we take in our simulations the standard Hasenfratz et al. values ai = 0.75, ^2 = 0.6, 
"3 = 0.3 Q- 



4. Renormalization factor Zb^^ 

The 4-quark operator renormalization needed for through one-loop can be written in 
terms of integrals that appear in the renormalization of bilinear operators. We thus calculate 
the renormalization factors for the quark bilinear operators O = qTq. The bilinear operator 
gets renormalized in the MS scheme according to 

Of^ifi) = (1 - wl)-'Z-'uoZr{fia)0'r^\l/a), (A33) 
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where Z^j renormalizes the domain- wall height. It is convenient to define the quark wave- 
function renormalization factor Z2 implicitly via the relation 



= (1 - w,)-'I^Z-JI\u,Z,fl\''^\ (A34) 



Using the Feynman rules presented in the previous sub-sections, we then have for the vertex 



correction to the bilinear operator in the MS*, NDR scheme 43 1 



Z2 lOvr 

with 



1 + ^ [Ar ln(/ia)2 + Ar(l - In vr^) + 5r - levr^Jr] , (A35) 



Ar = ^, i?^ = -^ + V^^^^ (A36) 

where h^i^) = A{V,A), 16{P,S), 0(r); = -1/2{V,A), 2{P,S), 0(T); and Ir is a finite 
lattice integral, 



^ /" J 5^ Tr [L{s)V,,{0, k) {i,{k, s)q{-k))T{q{k)^{-k, t))VU-k, 0)R{t)T^] 



X h^^{k)D^,{k)K,{k) - Ag'CpAr^^^JY^ ^ (A37) 



{k 

with the trace over Dirac spin and 



k 



"^'^ (A38) 



(2vr)4 

The last term in Eq. ( 1A37|) subtracts an IR (infrared) divergence from the integral. By 
chiral symmetry, the renormalization factors for the vector and axial-vector currents are 



equal; the renorma 



chiral symmetry 



ization factors for scalar and pseudoscalar currents are also equal by 



4l| . The Feynman diagram for the vertex correction is given in Fig. [31 



The renormalization factor matching the lattice calculation of Bk to the MS scheme can 



be written 



43| 



^"-^^"^ - (1 - " ^ ^ 

where Z^ is the renormalization factor for the operator O'^^^^, and Za renormalizes the axial 
current. It is useful to define in this way, since the tadpole and self-energy corrections 
cancel. The renormalization factor contains the running of the operator from the lattice 
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scale a ^ to the continuum scale fi. In the MS scheme with naive dimensional regularization 
(NDR), we obtain [431 



7MS,NDR/ 



a. 



MS,NDR 
Bk 



where 



z 



MS.NDR 
Bk 



with Isy defined in Eq. ( 1A37I) . 



(A40) 



(A41) 



Appendix B: Matching Scheme and Perturbative Running for Zbi. 



Although the functions used to convert the renormalization factor Zbj^ from the R 
scheme to the MS scheme are the same as those shown in the appendices of Ref. 
display them here for completeness. 



/MOM 



49|, we 



1. The QCD /3-function in the MS scheme 



In this work we calculate the value of the coupling constant (fi) at any scale using 
the four-loop (NNNLO) running formula of Ref. 77l |: 

d 



din 



a. 



where 



TT 



/3o 
Pi 

/S3 



1 / 2 

-fl02-^n, 
16 V 3 ^ ' ' 



1 /2857 5033 



64 



256 



149753 
6 



-Uf 



3564C3 



325 
'54 



-n, 



f ' 



1078361 6508 



+ 



/ 50065 6472 



V 162 



+ 



C3 n^f + 



162 
2 1093 ■ 
^ 729 ^ 



27 



(B2) 
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We implement this numerically by starting with the world average of the strong coupling 
constant at the Z-boson mass 



(m^) = 0.1176 ±0.0020, (B3) 

where the superscript indicates that this is determined in the region with five active quark 
fiavors. We then run ag below the bottom and charm quark thresholds imposing the match- 
ing conditions 

af) (mb) = af^ {nib) and af'^ (m,) = af^ (m,) , (B4) 
in order to determine al (/u) at any scale in the 3-fiavor theory. 



r: 

3 



2. Perturbative Running and Scheme Matching for Zbj^ 

We convert the renormalization factor Zbj^ between the scale- invariant, MS, and 
/MOM schemes using the one- loop renormalization group running formulae with Uf = 

Zt M = ^"Lme nf) (/X, nj) , (B5) 



where 



and 



1 (,, „ \ - ^MS / \-7o/2/3o 
heme l/^' ^f) " \N 



II) 

scheme 



^ ^ 4^ -^schen 



(B6) 



^nj j 17397 - 2070^/ + 104n2 

4i/MOM = 6(33-2n,f + ^ ' ' ^^^^ 



^(n,) ^ 13095 - 1626n^ + 8n} 

6(33-2r2/)' ■ ^ ' 

Appendix C: Non- perturbative Mixing Coefficients 

In this appendix we present results for the mixing coefficients between the operator Ok 
and other lattice operators of different chiralities. We compute them nonperturbatively 
using the method of Rome-Southampton as discussed in Sec. 1111 B II 

The renormalized operator that contributes to Bk in the continuum, which has a VV+AA 
chiral structure, receives contributions from several lattice operators: 

= Zvv+AA,i , (CI) 
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where i e {VV + AA, VV - AA, SS - PP, SS + PP, TT}. Because the operator mixings 
require two flips of chirality, the off-diagonal coefficients are suppressed by 0{{am 



res I I 



49 



52| . which is ~ 10"^ on the coarse lattice and even smaller on the fine lattice. We therefore 



expect the contributions to Bj^ from wrong-chirality lattice operators to be negligible. 

In practice, however, we find that the mixing coefficients are not of O {{arrij-csY) when we 
compute them using exceptional kinematics. This is because the choice of external momenta 
in the renormalization factor calculation leads to additional chiral symmetry breaking, as 
discussed in section [111 B 2[ Figures [T514T8] show the mixing coefficients as a function of (ap)^ 
for five valence quark masses on the ami/arrih = 0.007/0.05 coarse ensemble and in the 
chiral limit. At p ~ 2GeV, the mixing coefficients are still all quite small compared to Zbj^^- 
The largest is the mixing of Ok with the VV — AA operator, which is ~ 0.01. We observe 
coefficients of approximately the same size on the fine lattice, since this effect is not due to 
the lattice spacing or residual quark mass. 

Although the size of the mixing coefficients as computed with exceptional kinematics 
is not negligible, the results are contaminated by chiral symmetry breaking effects and 
are potentially unreliable. We therefore repeat the mixing coefficient calculation using non- 
exceptional kinematics. The results are shown for the coarse lattice in Figs. [T9] - [22] and for the 
fine lattice in Figs. [2314261 Although the mixing coefficients determined with non-exceptional 
kinematics have larger statistical errors, their values are smaller than when determined with 
exceptional kinematics. Furthermore, the new mixing coefficients are consistent with zero in 
the chiral limit. This confirms the hypothesis that the source of the large mixing coefficients 
is simply the choice of non-exceptional kinematics, and that the sizes of the true mixing 
coefficients are consistent with theoretical estimates. 

Although we find that the mixing coefficients are consistent with zero in the RI/MOM 
scheme using non-exceptional kinematics, we can still use this information to aid in our 
determination of Zbj^ using exceptional kinematics. This is because, ultimately, irrespective 
of the lattice scheme used to obtain the mixing coefficients, one must obtain the same mixing 
coefficients once the results are converted to the MS scheme. A vanishing contribution to 
BK from a particular operator in the RI/MOM scheme with non-exceptional kinematics 
implies a vanishing contribution in the MS scheme, since they are related multiplicatively. 
Generically, once an operator's contribution is zero in any scheme, its contribution is zero 
in all schemes that are multiplicatively related. Note, however, that once an operator's 
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am,/am^ = 0.007/0.05 

1 n 




0.5 1 1.5 2 2.5 3 

(ap)' 



FIG. 15: Mixing coefficient Zyyj^j^^Ayv-AA versus (ap)^ at several valence quark masses on the 
ami/arrih = 0.007/0.05 coarse ensemble. The stars indicate the value of the mixing coefficient in 
the limit that the valence quark mass and the light sea quark mass go to zero. 
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FIG. 18: Mixing coefficient Zvv+AA,TT versus (ap)^ at several valence quark masses on the 
ami/am^ = 0.007/0.05 coarse ensemble. 
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FIG. 24: Mixing coefficient Zvv+aa,ss~pp versus (ap)'^ at several valence quark masses on the 
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FIG. 25: Mixing coefficient Zyv+aa,ss+pp versus {ap)'^ at several valence quark masses on the 
ami/aiTih = 0.0062/0.031 fine ensemble, using non-exceptional kinematics. 
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FIG. 26: Mixing coefficient Zvv+aa,tt versus (ap)^ at several valence quark masses on the 
ami/aiTih = 0.0062/0.031 fine ensemble, using non-exceptional kinematics. 
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